1.資料介紹與下載

1.2資料下載:

  • 進入環保署環境資料開放平臺

  • 利用清單關鍵字收尋:「空氣品質監測小時值(一般污染物,每日更新)」

  • 選擇代碼為AQX_P_13

  • 下滑至歷史資料下載區塊

  • 勾選要下載的年-月檔案

  • 點選「下載勾選項目」,選擇csv檔(格式)

  • 下載後解壓縮得到csv 格式的檔案

  • 利用以下程式碼可將檔案匯入RStudio裡

  1. 設定工作區域
setwd("E:/rumiworkspace")
  1. 讀取資料
filename = "./data/G03_P13_2014_2022/空氣品質監測小時值資料(一般污染物,每日更新) (2014-01).csv"
df <- read.csv(filename,fileEncoding = "UTF-8-BOM")
  1. 檢查資料
head(df)
##   X.siteid. X.sitename. X.itemid. X.itemname. X.itemengname. X.itemunit.
## 1        11        士林         5    氮氧化物            NOx         ppb
## 2        11        士林         6    一氧化氮             NO         ppb
## 3        25        頭份        33  細懸浮微粒          PM2.5       μg/m3
## 4        25        頭份        14        溫度       AMB_TEMP           ℃
## 5        11        士林         7    二氧化氮            NO2         ppb
## 6         3        萬里        11        風向     WIND_DIREC     degrees
##   X.monitordate. X.monitorvalue00. X.monitorvalue01. X.monitorvalue02.
## 1     2014-01-01                16                15                29
## 2     2014-01-01               1.2                .8               3.3
## 3     2014-01-01                57                62                51
## 4     2014-01-01                14                13                13
## 5     2014-01-01                15                14                25
## 6     2014-01-01               212               210               213
##   X.monitorvalue03. X.monitorvalue04. X.monitorvalue05. X.monitorvalue06.
## 1                86                89                93               101
## 2                43                49                56                64
## 3                40                33                33                33
## 4                12                12                12                11
## 5                43                40                37                36
## 6               212               213               213               211
##   X.monitorvalue07. X.monitorvalue08. X.monitorvalue09. X.monitorvalue10.
## 1               119               129              20.0              14.0
## 2                82                82               5.7               3.3
## 3                38                39              34.0              31.0
## 4                11                14              17.0              20.0
## 5                37                47              15.0              11.0
## 6               209               192              69.0              44.0
##   X.monitorvalue11. X.monitorvalue12. X.monitorvalue13. X.monitorvalue14.
## 1                11                16                12                15
## 2               2.9               3.8               2.9               2.7
## 3                30                35                37                35
## 4                21                22                21                22
## 5               8.5                12               9.5                12
## 6                58                49                50                48
##   X.monitorvalue15. X.monitorvalue16. X.monitorvalue17. X.monitorvalue18.
## 1                16                14                17                14
## 2               2.5               2.3               1.9               1.9
## 3                38                45                51                59
## 4                22                21                20                19
## 5                14                12                15                12
## 6                57                77                66               200
##   X.monitorvalue19. X.monitorvalue20. X.monitorvalue21. X.monitorvalue22.
## 1                13                11                11               8.6
## 2               1.7               1.9               1.9               1.6
## 3                65                62                56              51.0
## 4                18                17                17              17.0
## 5                11               8.7                 9               6.9
## 6               208               205               206             200.0
##   X.monitorvalue23.
## 1               6.8
## 2               1.7
## 3              47.0
## 4              17.0
## 5               5.1
## 6             207.0

1.3變項介紹:

•SiteId(測站代碼):

每個測站都有一個自己的代碼,代表該測站的識別(1~72、75、77、78、80、83、84、85)

•SiteName(測站名稱):

是測站的名稱(基隆、汐止、萬里、…、大城)

•ItemId(測項代碼):

每個測項都有一個自己的代碼,代表該測項的識別(1~9、11、14、21~23、31~33、36、38、143~144)

•ItemName(測項名稱):

是測項的名稱(二氧化硫、一氧化碳、臭氧、…、小時風向值)

•ItemEngName(測項英文名稱):

是測項的英文名稱(SO2、CO、O3、…、WD_HR)

•ItemUnit(測項單位):

測項的計量單位,如 μg/m³、ppm等。

Item表格整理+說明✔

ItemId ItemName ItemEngName ItemUnit 測項說明
1 二氧化硫 SO2 ppb 空氣中的一種污染物,主要來自燃煤、石油等化石燃料的燃燒排放。
2 一氧化碳 CO ppm 空氣中的一種有毒氣體,主要來自車輛排放和燃燒過程中不完全燃燒。
3 臭氧 O3 ppb 空氣中的一種臭氣,由氮氧化物和揮發性有機物在陽光照射下產生,對人體和植物有害。
4 懸浮微粒 PM10 μg/m3 空氣中的固體和液體顆粒物,可造成呼吸道、心血管等健康問題。
5 氮氧化物 NOx ppb 包括一氧化氮和二氧化氮,主要來自燃煤、汽車等排放源。
6 一氧化氮 NO ppb 一種有毒氣體,主要來自汽車排放和燃燒過程中高溫氮氧化物的生成。
7 二氧化氮 NO2 ppb 一種有毒氣體,主要來自汽車、工廠等排放源,可造成呼吸道問題。
8 總碳氫化物 THC ppm 包括甲烷、乙烷、丙烷等氣體,以及揮發性有機化合物,主要來自燃料燃燒和工業過程。
9 非甲烷碳氫化合物 NMHC ppm 不包括甲烷的碳氫化合物,主要來自石油和天然氣開采、運輸等過程。
11 風向 WIND_DIREC degrees 風的方向。
14 溫度 AMB_TEMP 空氣的溫度,通常用攝氏度或華氏度表示。
21 酸雨 PH_RAIN pH 指大氣中酸性物質與雨水混合而成的雨水,對生態環境和建築物有害。
22 導電度 RAIN_COND μmho/cm 水體中離子的濃度和移動能力的度量,通常用西門子/厘米表示。
23 雨量 RAINFALL mm 一定時間內的雨水總量,通常用毫米表示。
31 甲烷 CH4 ppm 一種無色、無味、易燃的氣體,是溫室氣體之一。
32 降雨強度 RAIN_INT 降雨時單位時間內雨量的多寡,常用毫米每小時(mm/hr)來表示。
33 細懸浮微粒 PM2.5 μg/m3 空氣中直徑小於2.5微米的顆粒物,會對人體健康造成危害,如引發呼吸道疾病。
36 二氧化碳 CO2 ppm 一種無色、無味、不易燃的氣體,是主要的溫室氣體之一,對全球氣候變化有重要影響。
38 相對濕度 RH % 空氣中水蒸氣的含量與空氣中最大可能的含量之比,通常以百分比表示。
143 小時風速值 WS_HR m/sec 一小時內風速的最大值。
144 小時風向值 WD_HR degrees 一小時內風向的平均值。

•MonitorDate(監測日期):

測量數據的日期(20xx/xx/xx)

•MonitorValue0(0時數值)~MonitorValue23(23時數值):

測量數據的具體數值,代表在該時的空氣品質監測值,單位為測項的計量

1.4變項之間相關性介紹:

  1. SiteId 和 SiteName 是相關的,表示同一個測站

    在分析數據時,可以根據 SiteName 篩選出特定城市或地區的監測數據,或者根據 SiteId 查詢該測站的歷史數據。

  2. MonitorDate 和 MonitorValue0~MonitorValue23 是相關的,表示同一天的不同時間點的監測數據

    可以分析數據的時間趨勢和變化情況。

  3. ItemId 和 ItemName/ItemEngName:每個 ItemId 都對應一個 ItemName 和一個 ItemEngName,表示測量的物質名稱和英文名稱,因此這三個變項是相關的。

example✔

在分析數據時,可以根據 ItemName/ItemEngName 篩選出特定物質的監測數據,或者根據 ItemId 查詢該測項的歷史數據。

setwd("E:/rumiworkspace")
filename = "./data/datatest/空氣品質監測小時值資料(一般污染物,每日更新)_test.csv"
df = read.csv(filename,fileEncoding = "UTF-8-BOM")
datatable(head(df, 100), options = list(scrollx = TRUE, scrollCollapse = TRUE))

2.資料整理

目標 利用 R 的 dplyr 的套件中的指令,來整理資料,並且轉置資料

2.1 了解並載入dplyr套件

#載入dplyr套件
library(dplyr)

補充:淺談dplyr基礎功能

  • select():從數據中選擇列

  • filter():數據行的子集 (篩選符合條件的觀察值)

  • group_by():匯總數據 (根據資料內的類別變數對另一個變數進行統計)

  • summarise():匯總數據(計算匯總統計信息)

  • arrange():排序數據 (變數遞增或是遞減排序資料)

  • mutate():創建新變量

示範:

#自訂新變量[newcode]
setwd("E:/rumiworkspace")
airpollution<-read.csv('./data/datatest/空氣品質監測小時值資料(一般污染物,每日更新)_test.csv'
,fileEncoding = "UTF-8-BOM")
airpollution.part<-mutate(airpollution, newcode=round(as.numeric(averageday)))
datatable(head(airpollution.part, 100), options = list(scrollx = TRUE, scrollCollapse =
TRUE))
#利用filter可以抓取指定測項和其規範數值
setwd("E:/rumiworkspace")
airpollution<-read.csv('./data/datatest/空氣品質監測小時值資料(一般污染物,每日更新)_test.csv'
,fileEncoding = "UTF-8-BOM")
airpollution.filter<-filter(airpollution, itemengname == 'NOx',as.numeric(averageday) > 
100)
airpollution.filter
##    itemname itemengname itemunit monitordate averageday
## 1  氮氧化物         NOx      ppb    2014/1/1     135.63
## 2  氮氧化物         NOx      ppb    2014/1/2     110.00
## 3  氮氧化物         NOx      ppb    2014/1/2     194.73
## 4  氮氧化物         NOx      ppb    2014/1/2     100.98
## 5  氮氧化物         NOx      ppb    2014/1/2     129.25
## 6  氮氧化物         NOx      ppb    2014/1/3     174.67
## 7  氮氧化物         NOx      ppb    2014/1/3     106.13
## 8  氮氧化物         NOx      ppb    2014/1/4     142.96
## 9  氮氧化物         NOx      ppb    2014/1/5     113.38
## 10 氮氧化物         NOx      ppb    2014/1/6     124.17
## 11 氮氧化物         NOx      ppb    2014/1/7     130.46
## 12 氮氧化物         NOx      ppb    2014/1/7     102.29
## 13 氮氧化物         NOx      ppb    2014/1/7     193.88
## 14 氮氧化物         NOx      ppb    2014/1/8     110.54
## 15 氮氧化物         NOx      ppb    2014/1/8     187.17
## 16 氮氧化物         NOx      ppb    2014/1/9     124.88
## 17 氮氧化物         NOx      ppb   2014/1/11     119.13
## 18 氮氧化物         NOx      ppb   2014/1/12     124.88

2.2 資料察看

以了解我們資料是否有地方可以修改。

# 檔案名稱
setwd("E:/rumiworkspace")
filename = "./data/G03_P13_2014_2022/空氣品質監測小時值資料(一般污染物,每日更新) (2014-01).csv"
# 讀入資料
df = read.csv(filename, fileEncoding="UTF-8-BOM")
str(df)
## 'data.frame':    27756 obs. of  31 variables:
##  $ X.siteid.        : int  11 11 25 25 11 3 42 25 71 11 ...
##  $ X.sitename.      : chr  "士林" "士林" "頭份" "頭份" ...
##  $ X.itemid.        : int  5 6 33 14 7 11 4 11 10 10 ...
##  $ X.itemname.      : chr  "氮氧化物" "一氧化氮" "細懸浮微粒" "溫度" ...
##  $ X.itemengname.   : chr  "NOx" "NO" "PM2.5" "AMB_TEMP" ...
##  $ X.itemunit.      : chr  "ppb" "ppb" "μg/m3" "℃" ...
##  $ X.monitordate.   : chr  "2014-01-01" "2014-01-01" "2014-01-01" "2014-01-01" ...
##  $ X.monitorvalue00.: chr  "16" "1.2" "57" "14" ...
##  $ X.monitorvalue01.: chr  "15" ".8" "62" "13" ...
##  $ X.monitorvalue02.: chr  "29" "3.3" "51" "13" ...
##  $ X.monitorvalue03.: chr  "86" "43" "40" "12" ...
##  $ X.monitorvalue04.: chr  "89" "49" "33" "12" ...
##  $ X.monitorvalue05.: chr  "93" "56" "33" "12" ...
##  $ X.monitorvalue06.: chr  "101" "64" "33" "11" ...
##  $ X.monitorvalue07.: chr  "119" "82" "38" "11" ...
##  $ X.monitorvalue08.: chr  "129" "82" "39" "14" ...
##  $ X.monitorvalue09.: chr  "20.0" "5.7" "34.0" "17.0" ...
##  $ X.monitorvalue10.: chr  "14.0" "3.3" "31.0" "20.0" ...
##  $ X.monitorvalue11.: chr  "11" "2.9" "30" "21" ...
##  $ X.monitorvalue12.: chr  "16" "3.8" "35" "22" ...
##  $ X.monitorvalue13.: chr  "12" "2.9" "37" "21" ...
##  $ X.monitorvalue14.: chr  "15" "2.7" "35" "22" ...
##  $ X.monitorvalue15.: chr  "16" "2.5" "38" "22" ...
##  $ X.monitorvalue16.: chr  "14" "2.3" "45" "21" ...
##  $ X.monitorvalue17.: chr  "17" "1.9" "51" "20" ...
##  $ X.monitorvalue18.: chr  "14" "1.9" "59" "19" ...
##  $ X.monitorvalue19.: chr  "13" "1.7" "65" "18" ...
##  $ X.monitorvalue20.: chr  "11" "1.9" "62" "17" ...
##  $ X.monitorvalue21.: chr  "11" "1.9" "56" "17" ...
##  $ X.monitorvalue22.: chr  "8.6" "1.6" "51.0" "17.0" ...
##  $ X.monitorvalue23.: chr  "6.8" "1.7" "47.0" "17.0" ...
head(df)
##   X.siteid. X.sitename. X.itemid. X.itemname. X.itemengname. X.itemunit.
## 1        11        士林         5    氮氧化物            NOx         ppb
## 2        11        士林         6    一氧化氮             NO         ppb
## 3        25        頭份        33  細懸浮微粒          PM2.5       μg/m3
## 4        25        頭份        14        溫度       AMB_TEMP           ℃
## 5        11        士林         7    二氧化氮            NO2         ppb
## 6         3        萬里        11        風向     WIND_DIREC     degrees
##   X.monitordate. X.monitorvalue00. X.monitorvalue01. X.monitorvalue02.
## 1     2014-01-01                16                15                29
## 2     2014-01-01               1.2                .8               3.3
## 3     2014-01-01                57                62                51
## 4     2014-01-01                14                13                13
## 5     2014-01-01                15                14                25
## 6     2014-01-01               212               210               213
##   X.monitorvalue03. X.monitorvalue04. X.monitorvalue05. X.monitorvalue06.
## 1                86                89                93               101
## 2                43                49                56                64
## 3                40                33                33                33
## 4                12                12                12                11
## 5                43                40                37                36
## 6               212               213               213               211
##   X.monitorvalue07. X.monitorvalue08. X.monitorvalue09. X.monitorvalue10.
## 1               119               129              20.0              14.0
## 2                82                82               5.7               3.3
## 3                38                39              34.0              31.0
## 4                11                14              17.0              20.0
## 5                37                47              15.0              11.0
## 6               209               192              69.0              44.0
##   X.monitorvalue11. X.monitorvalue12. X.monitorvalue13. X.monitorvalue14.
## 1                11                16                12                15
## 2               2.9               3.8               2.9               2.7
## 3                30                35                37                35
## 4                21                22                21                22
## 5               8.5                12               9.5                12
## 6                58                49                50                48
##   X.monitorvalue15. X.monitorvalue16. X.monitorvalue17. X.monitorvalue18.
## 1                16                14                17                14
## 2               2.5               2.3               1.9               1.9
## 3                38                45                51                59
## 4                22                21                20                19
## 5                14                12                15                12
## 6                57                77                66               200
##   X.monitorvalue19. X.monitorvalue20. X.monitorvalue21. X.monitorvalue22.
## 1                13                11                11               8.6
## 2               1.7               1.9               1.9               1.6
## 3                65                62                56              51.0
## 4                18                17                17              17.0
## 5                11               8.7                 9               6.9
## 6               208               205               206             200.0
##   X.monitorvalue23.
## 1               6.8
## 2               1.7
## 3              47.0
## 4              17.0
## 5               5.1
## 6             207.0

從以上兩表發現有幾處是需要修改的。

  1. 欄位名稱需修改
原名稱 英文名稱 中文名稱
X.siteid. siteid 測站代碼
X.sitename. sitename 測站名稱(V)
X.itemid. itemid 測項代碼
X.itemname. itemname 測項名稱
X.itemengname. itemengname 測項英文名稱(V)
X.itemunit. itemunit 測項單位
X.monitordate. monitordate 測量日期(V)
X.monitorvalue00. ~ 23 h00~h23 測量值(V)
  1. 欄位資料類型需要修改

    • monitordate須為日期資料類型

    • 測量值需為數值資料類型

2.3 測站資訊

查看測站相關變數。

有下列資訊:

  • 測站代碼(siteid)

  • 測站名稱(sitename)

# 測站資訊檔
stateInfo = data.frame(df$X.siteid., df$X.sitename.)
colnames(stateInfo) = c("測站代碼", "測站名稱")
stateInfo = unique(stateInfo)
str(stateInfo)
## 'data.frame':    76 obs. of  2 variables:
##  $ 測站代碼: int  11 25 3 42 71 35 56 72 12 4 ...
##  $ 測站名稱: chr  "士林" "頭份" "萬里" "嘉義" ...

2.4 測項資訊

查看測項相關變數。

有下列資訊:

  • 測項代碼(itemid)

  • 測項中文名稱(itemname)

  • 測項英文名稱(itemengname)

  • 測項單位(itemunit)

# 測項資訊檔
itemInfo = data.frame(df$X.itemid., df$X.itemname., 
                      df$X.itemengname., df$X.itemunit.)
colnames(itemInfo) = c("測項代碼", "測項中文名稱", "測項英文名稱", "測項單位")
itemInfo = unique(itemInfo)
str(itemInfo)
## 'data.frame':    12 obs. of  4 variables:
##  $ 測項代碼    : int  5 6 33 14 7 11 4 10 3 38 ...
##  $ 測項中文名稱: chr  "氮氧化物" "一氧化氮" "細懸浮微粒" "溫度" ...
##  $ 測項英文名稱: chr  "NOx" "NO" "PM2.5" "AMB_TEMP" ...
##  $ 測項單位    : chr  "ppb" "ppb" "μg/m3" "℃" ...

2.5 測量資料(欄位名稱)更改

以便於我們識別欄位。

# 更改欄位名稱
dataInfo = dplyr::select(
  df, X.sitename., X.itemengname., 
  X.monitordate., starts_with("X.monitorvalue")
)

colnames(dataInfo) = c("測站名稱", "測項英文名稱", "測量日期", 
                     paste("h0",seq(0,9),sep=""),
                     paste("h",seq(10,23),sep=""))
str(dataInfo)
## 'data.frame':    27756 obs. of  27 variables:
##  $ 測站名稱    : chr  "士林" "士林" "頭份" "頭份" ...
##  $ 測項英文名稱: chr  "NOx" "NO" "PM2.5" "AMB_TEMP" ...
##  $ 測量日期    : chr  "2014-01-01" "2014-01-01" "2014-01-01" "2014-01-01" ...
##  $ h00         : chr  "16" "1.2" "57" "14" ...
##  $ h01         : chr  "15" ".8" "62" "13" ...
##  $ h02         : chr  "29" "3.3" "51" "13" ...
##  $ h03         : chr  "86" "43" "40" "12" ...
##  $ h04         : chr  "89" "49" "33" "12" ...
##  $ h05         : chr  "93" "56" "33" "12" ...
##  $ h06         : chr  "101" "64" "33" "11" ...
##  $ h07         : chr  "119" "82" "38" "11" ...
##  $ h08         : chr  "129" "82" "39" "14" ...
##  $ h09         : chr  "20.0" "5.7" "34.0" "17.0" ...
##  $ h10         : chr  "14.0" "3.3" "31.0" "20.0" ...
##  $ h11         : chr  "11" "2.9" "30" "21" ...
##  $ h12         : chr  "16" "3.8" "35" "22" ...
##  $ h13         : chr  "12" "2.9" "37" "21" ...
##  $ h14         : chr  "15" "2.7" "35" "22" ...
##  $ h15         : chr  "16" "2.5" "38" "22" ...
##  $ h16         : chr  "14" "2.3" "45" "21" ...
##  $ h17         : chr  "17" "1.9" "51" "20" ...
##  $ h18         : chr  "14" "1.9" "59" "19" ...
##  $ h19         : chr  "13" "1.7" "65" "18" ...
##  $ h20         : chr  "11" "1.9" "62" "17" ...
##  $ h21         : chr  "11" "1.9" "56" "17" ...
##  $ h22         : chr  "8.6" "1.6" "51.0" "17.0" ...
##  $ h23         : chr  "6.8" "1.7" "47.0" "17.0" ...

2.6 轉置資料(寬表與長表)

何謂轉置?

在線性代數的矩陣操作裡,依據左上至右下的對角線來進行列資料與欄資料的對換稱為「轉置」(Transpose);在表格式資料(試算表、關聯表格或資料框)中也有相似的操作,由於表格式資料相對於矩陣,多出了變數的欄位名稱與觀測值的列索引,轉置前後的外觀又被稱作寬表(Wide Format)與長表(Long Format)。

為何要轉置資料呢?

資料的外觀與函式並無法呼應。這時資料轉置的需求就相當地明確

我們有兩種欄位類型:

  • 長表(long):測站名稱,測項名稱,測量日期,測量小時,測量值
# 轉製成長版
varying = c(paste("h0",seq(0,9),sep=""),
          paste("h",seq(10,23),sep=""))
df.long = reshape(dataInfo, varying=varying, 
              v.names="測量值",
              timevar="測量小時",
              #times=varying, #seq_along(varying[[1]]),
              direction="long") 
str(df.long)
## 'data.frame':    666144 obs. of  6 variables:
##  $ 測站名稱    : chr  "士林" "士林" "頭份" "頭份" ...
##  $ 測項英文名稱: chr  "NOx" "NO" "PM2.5" "AMB_TEMP" ...
##  $ 測量日期    : chr  "2014-01-01" "2014-01-01" "2014-01-01" "2014-01-01" ...
##  $ 測量小時    : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ 測量值      : chr  "16" "1.2" "57" "14" ...
##  $ id          : int  1 2 3 4 5 6 7 8 9 10 ...
##  - attr(*, "reshapeLong")=List of 4
##   ..$ varying:List of 1
##   .. ..$ 測量值: chr [1:24] "h00" "h01" "h02" "h03" ...
##   .. ..- attr(*, "v.names")= chr "測量值"
##   .. ..- attr(*, "times")= int [1:24] 1 2 3 4 5 6 7 8 9 10 ...
##   ..$ v.names: chr "測量值"
##   ..$ idvar  : chr "id"
##   ..$ timevar: chr "測量小時"
df.long$id = NULL
attr(df.long, "reshapeLong") = NULL
str(df.long)
## 'data.frame':    666144 obs. of  5 variables:
##  $ 測站名稱    : chr  "士林" "士林" "頭份" "頭份" ...
##  $ 測項英文名稱: chr  "NOx" "NO" "PM2.5" "AMB_TEMP" ...
##  $ 測量日期    : chr  "2014-01-01" "2014-01-01" "2014-01-01" "2014-01-01" ...
##  $ 測量小時    : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ 測量值      : chr  "16" "1.2" "57" "14" ...
head(df.long, n=30)
##      測站名稱 測項英文名稱   測量日期 測量小時 測量值
## 1.1      士林          NOx 2014-01-01        1     16
## 2.1      士林           NO 2014-01-01        1    1.2
## 3.1      頭份        PM2.5 2014-01-01        1     57
## 4.1      頭份     AMB_TEMP 2014-01-01        1     14
## 5.1      士林          NO2 2014-01-01        1     15
## 6.1      萬里   WIND_DIREC 2014-01-01        1    212
## 7.1      嘉義         PM10 2014-01-01        1     91
## 8.1      頭份   WIND_DIREC 2014-01-01        1     71
## 9.1      復興   WIND_SPEED 2014-01-01        1     .8
## 10.1     士林   WIND_SPEED 2014-01-01        1    1.2
## 11.1     二林          NOx 2014-01-01        1     23
## 12.1     士林         PM10 2014-01-01        1     56
## 13.1     復興   WIND_DIREC 2014-01-01        1    113
## 14.1     頭份   WIND_SPEED 2014-01-01        1     .9
## 15.1     士林   WIND_DIREC 2014-01-01        1    141
## 16.1     萬里     AMB_TEMP 2014-01-01        1     13
## 17.1     頭份          NO2 2014-01-01        1     24
## 18.1     復興     AMB_TEMP 2014-01-01        1     17
## 19.1     士林     AMB_TEMP 2014-01-01        1     16
## 20.1     頭份           NO 2014-01-01        1      2
## 21.1     士林           O3 2014-01-01        1     45
## 22.1     復興        PM2.5 2014-01-01        1     83
## 23.1     復興           RH 2014-01-01        1     74
## 24.1     士林        PM2.5 2014-01-01        1     28
## 25.1     前金          NO2 2014-01-01        1     30
## 26.1     頭份          NOx 2014-01-01        1     26
## 27.1     萬里        PM2.5 2014-01-01        1     30
## 28.1     士林           RH 2014-01-01        1     67
## 29.1     埔里          SO2 2014-01-01        1    2.5
## 30.1     二林         PM10 2014-01-01        1     78
  • 資料類型轉換

在寬版轉置前,需要進行2.2 資料察看的第2.點->「資料轉換」

首先,通常在資料整理時,必須參考資料來源單位所提供的說明檔。

在說明檔中有下列針對普通測站的說明:

  • # 表示儀器檢核為無效值

  • * 表示程式檢核為無效值

  • x 表示人工檢核為無效值

  • NR 表示無降雨

  • 空白 表示缺值

  • 風向資料888代表無風,999則代表儀器故障。

我可以利用grepl()來查看是否有特殊字串:

q1 = grepl("#", df.long$測量值, fixed=T)
sum(q1)
## [1] 0
q1 = grepl("*", df.long$測量值, fixed=T)
sum(q1)
## [1] 0
q1 = grepl("x", df.long$測量值, fixed=T)
sum(q1)
## [1] 6720
q1 = grepl("NR", df.long$測量值, fixed=T)
sum(q1)
## [1] 0
q1 = grepl("888", df.long$測量值, fixed=T)
sum(q1)
## [1] 0
q1 = grepl("999", df.long$測量值, fixed=T)
sum(q1)
## [1] 0
sum(nchar(df.long$測量值)==0)
## [1] 0

沒有異樣,可以進行轉換:

# 資料轉換
df.long[, "測量值old"] = df.long[, "測量值"]
df.long[, "測量值"] = suppressWarnings(as.numeric(df.long[, "測量值"]))
summary(df.long[, -c(1,2,3)])
##     測量小時         測量值        測量值old        
##  Min.   : 1.00   Min.   : -1.00   Length:666144     
##  1st Qu.: 6.75   1st Qu.:  3.20   Class :character  
##  Median :12.50   Median : 18.00   Mode  :character  
##  Mean   :12.50   Mean   : 36.94                     
##  3rd Qu.:18.25   3rd Qu.: 50.00                     
##  Max.   :24.00   Max.   :968.00                     
##                  NA's   :9671
# 日期資料轉換
df.long[, "測量日期"] = as.Date(df.long[, "測量日期"])
str(df.long)
## 'data.frame':    666144 obs. of  6 variables:
##  $ 測站名稱    : chr  "士林" "士林" "頭份" "頭份" ...
##  $ 測項英文名稱: chr  "NOx" "NO" "PM2.5" "AMB_TEMP" ...
##  $ 測量日期    : Date, format: "2014-01-01" "2014-01-01" ...
##  $ 測量小時    : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ 測量值      : num  16 1.2 57 14 15 212 91 71 0.8 1.2 ...
##  $ 測量值old   : chr  "16" "1.2" "57" "14" ...
  • 寬表(widw):測站名稱,測量日期,測量小時,每個測量項目一欄
# 轉製成寬版
dfTemp = df.long
dfTemp[, "測量值old"] = NULL

idvar = c("測站名稱", "測量日期", "測量小時")
df.wide = reshape(dfTemp, idvar=idvar, 
              timevar="測項英文名稱",
              v.names="測量值",
              direction="wide") 
str(df.wide)
## 'data.frame':    56448 obs. of  15 variables:
##  $ 測站名稱         : chr  "士林" "頭份" "萬里" "嘉義" ...
##  $ 測量日期         : Date, format: "2014-01-01" "2014-01-01" ...
##  $ 測量小時         : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ 測量值.NOx       : num  16 26 5 37 75 23 32 17 46 20 ...
##  $ 測量值.NO        : num  1.2 2 0.9 3.1 29 1.4 1.3 1 5.8 2.5 ...
##  $ 測量值.PM2.5     : num  28 57 30 65 83 46 85 69 59 24 ...
##  $ 測量值.AMB_TEMP  : num  16 14 13 15 17 13 16 14 17 15 ...
##  $ 測量值.NO2       : num  15 24 4.1 34 46 21 30 16 41 17 ...
##  $ 測量值.WIND_DIREC: num  141 71 212 61 113 38 34 312 128 118 ...
##  $ 測量值.PM10      : num  56 76 38 91 134 78 120 79 76 40 ...
##  $ 測量值.WIND_SPEED: num  1.2 0.9 4.7 1.2 0.8 1.8 1.7 0.9 1.5 0.9 ...
##  $ 測量值.O3        : num  45 7.9 36 7.7 1.5 20 14 17 13 22 ...
##  $ 測量值.RH        : num  67 80 86 80 74 81 78 80 70 76 ...
##  $ 測量值.SO2       : num  6.3 2.8 1.7 5.1 7.6 4.3 4.9 2.5 3.1 2 ...
##  $ 測量值.CO        : num  0.44 0.53 0.25 0.88 1.48 0.44 0.77 0.65 0.87 0.53 ...
##  - attr(*, "reshapeWide")=List of 5
##   ..$ v.names: chr "測量值"
##   ..$ timevar: chr "測項英文名稱"
##   ..$ idvar  : chr [1:3] "測站名稱" "測量日期" "測量小時"
##   ..$ times  : chr [1:12] "NOx" "NO" "PM2.5" "AMB_TEMP" ...
##   ..$ varying: chr [1, 1:12] "測量值.NOx" "測量值.NO" "測量值.PM2.5" "測量值.AMB_TEMP" ...
# 整理欄位名稱
colnames(df.wide)[-c(1,2,3)] = substr(colnames(df.wide)[-c(1,2,3)], 
                                  5, nchar(colnames(df.wide)[-c(1,2,3)]))
attr(df.wide, "reshapeWide") = NULL
str(df.wide)
## 'data.frame':    56448 obs. of  15 variables:
##  $ 測站名稱  : chr  "士林" "頭份" "萬里" "嘉義" ...
##  $ 測量日期  : Date, format: "2014-01-01" "2014-01-01" ...
##  $ 測量小時  : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ NOx       : num  16 26 5 37 75 23 32 17 46 20 ...
##  $ NO        : num  1.2 2 0.9 3.1 29 1.4 1.3 1 5.8 2.5 ...
##  $ PM2.5     : num  28 57 30 65 83 46 85 69 59 24 ...
##  $ AMB_TEMP  : num  16 14 13 15 17 13 16 14 17 15 ...
##  $ NO2       : num  15 24 4.1 34 46 21 30 16 41 17 ...
##  $ WIND_DIREC: num  141 71 212 61 113 38 34 312 128 118 ...
##  $ PM10      : num  56 76 38 91 134 78 120 79 76 40 ...
##  $ WIND_SPEED: num  1.2 0.9 4.7 1.2 0.8 1.8 1.7 0.9 1.5 0.9 ...
##  $ O3        : num  45 7.9 36 7.7 1.5 20 14 17 13 22 ...
##  $ RH        : num  67 80 86 80 74 81 78 80 70 76 ...
##  $ SO2       : num  6.3 2.8 1.7 5.1 7.6 4.3 4.9 2.5 3.1 2 ...
##  $ CO        : num  0.44 0.53 0.25 0.88 1.48 0.44 0.77 0.65 0.87 0.53 ...
head(df.wide)
##      測站名稱   測量日期 測量小時 NOx   NO PM2.5 AMB_TEMP  NO2 WIND_DIREC PM10
## 1.1      士林 2014-01-01        1  16  1.2    28       16 15.0        141   56
## 3.1      頭份 2014-01-01        1  26  2.0    57       14 24.0         71   76
## 6.1      萬里 2014-01-01        1   5  0.9    30       13  4.1        212   38
## 7.1      嘉義 2014-01-01        1  37  3.1    65       15 34.0         61   91
## 9.1      復興 2014-01-01        1  75 29.0    83       17 46.0        113  134
## 11.1     二林 2014-01-01        1  23  1.4    46       13 21.0         38   78
##      WIND_SPEED   O3 RH SO2   CO
## 1.1         1.2 45.0 67 6.3 0.44
## 3.1         0.9  7.9 80 2.8 0.53
## 6.1         4.7 36.0 86 1.7 0.25
## 7.1         1.2  7.7 80 5.1 0.88
## 9.1         0.8  1.5 74 7.6 1.48
## 11.1        1.8 20.0 81 4.3 0.44

2.7 儲存檔案

# 輸出檔案名稱
outDir = "./data/abc"
p1 = gregexpr(".csv", filename)

# 利用檔名製作輸出檔名
outName = paste(outDir, "/L",
                substr(filename, p1[[1]][1]-8, p1[[1]][1]-2),
                ".rds", sep="")

saveRDS(df.long, outName)


# 輸出檔案名稱
outDir = "./data/abc"
p1 = gregexpr(".csv", filename)

# 利用檔名製作輸出檔名
outName = paste(outDir, "/W",
                substr(filename, p1[[1]][1]-8, p1[[1]][1]-2),
                ".rds", sep="")

saveRDS(df.wide, outName)

這樣檔案就輸出至abc資料夾裡了:

2.8 寫成函數(並大量作業)

由於有眾多檔案,我們不可能一個一個整理輸出,此時我們將上述程式碼整理成函數,使該函數可以大量作業,無須我們擠一個一個慢慢來。

程式碼如下:

pre1 = function(filename, outDir){
    
  # 載入套件
  library(dplyr)
  # 讀入資料
  df = read.csv(filename, fileEncoding="UTF-8-BOM")
  if (sum(nchar(df$X.monitordate.)!=10)) df$X.monitordate. = substr(df$X.monitordate., 1, 10)
  # 防止重複資料
  df = unique(df)
  # 測站資訊檔
  stateInfo = data.frame(df$X.siteid., df$X.sitename.)
  colnames(stateInfo) = c("測站代碼", "測站名稱")
  stateInfo = unique(stateInfo)
  # 測項資訊檔
  itemInfo = data.frame(df$X.itemid., df$X.itemname., 
                        df$X.itemengname., df$X.itemunit.)
  colnames(itemInfo) = c("測項代碼", "測項中文名稱", "測項英文名稱", "測項單位")
  itemInfo = unique(itemInfo)
  # 更改欄位名稱
  dataInfo = dplyr::select(
    df, X.sitename., X.itemengname., 
    X.monitordate., starts_with("X.monitorvalue")
  )
  colnames(dataInfo) = c("測站名稱", "測項英文名稱", "測量日期", 
                       paste("h0",seq(0,9),sep=""),
                       paste("h",seq(10,23),sep=""))
  # 轉製成長版
  varying = c(paste("h0",seq(0,9),sep=""),
            paste("h",seq(10,23),sep=""))
  df.long = reshape(dataInfo, varying=varying, 
                v.names="測量值",
                timevar="測量小時",
                direction="long") 
  df.long$id = NULL
  attr(df.long, "reshapeLong") = NULL
  df.long[, "測量值old"] = df.long[, "測量值"]
  df.long[, "測量值"] = suppressWarnings(as.numeric(df.long[, "測量值"]))
  # 日期資料轉換
  df.long[, "測量日期"] = as.Date(df.long[, "測量日期"])
  # 排序
  # 轉製成寬版
  dfTemp = df.long
  dfTemp[, "測量值old"] = NULL
  idvar = c("測站名稱", "測量日期", "測量小時")
  df.wide = reshape(dfTemp, idvar=idvar, 
                timevar="測項英文名稱",
                v.names="測量值",
                direction="wide")
  # 整理欄位名稱
  colnames(df.wide)[-c(1,2,3)] = substr(colnames(df.wide)[-c(1,2,3)], 
                                    5, nchar(colnames(df.wide)[-c(1,2,3)]))
  attr(df.wide, "reshapeWide") = NULL
  #===長板===
  # 利用檔名製作輸出檔名
  p1 = gregexpr(".csv", filename)
  outName = paste(outDir, "/L",
                  substr(filename, p1[[1]][1]-8, p1[[1]][1]-2),
                  ".rds", sep="")
  
  saveRDS(df.long, outName)
  
  #===寬板===
  # 輸出檔案名稱
  p1 = gregexpr(".csv", filename)
  
  # 利用檔名製作輸出檔名
  outName = paste(outDir, "/W",
                  substr(filename, p1[[1]][1]-8, p1[[1]][1]-2),
                  ".rds", sep="")
  
  saveRDS(df.wide, outName)

  # 輸出
  return(list(stateInfo=stateInfo, itemInfo=itemInfo))
}

測試 pre1() 函數:

filename = "./data/G03_P13_2014_2022/空氣品質監測小時值資料(一般污染物,每日更新) (2014-01).csv"

outDir="./data/S_0165_P13_G03"
q1 = pre1(filename, outDir=outDir)
head(q1$stateInfo)
##    測站代碼 測站名稱
## 1        11     士林
## 3        25     頭份
## 6         3     萬里
## 7        42     嘉義
## 9        71     復興
## 11       35     二林
q1$itemInfo
##    測項代碼 測項中文名稱 測項英文名稱 測項單位
## 1         5     氮氧化物          NOx      ppb
## 2         6     一氧化氮           NO      ppb
## 3        33   細懸浮微粒        PM2.5    μg/m3
## 4        14         溫度     AMB_TEMP        ℃
## 5         7     二氧化氮          NO2      ppb
## 6        11         風向   WIND_DIREC  degrees
## 7         4     懸浮微粒         PM10    μg/m3
## 9        10         風速   WIND_SPEED    m/sec
## 21        3         臭氧           O3      ppb
## 23       38     相對濕度           RH  percent
## 29        1     二氧化硫          SO2      ppb
## 31        2     一氧化碳           CO      ppm
  • 利用list.dirs()查詢出某一目錄下的全部目錄名稱。

    • 查出來的第一個目錄是本身。
  • 利用list.files()查詢出某一目錄下的全部檔案名稱。

接著,就可以進行2014-2022的大量轉置作業~

補充:由於共9年,每年12個月,所以要輸入的數字為[1:108] !

# 檔案所在目錄
target = "./data/G03_P13_2014_2022/"
stateTemp = NULL
itemTemp = NULL
outDir = "./data/S_0165_P13_G03"

  q2 = list.files(target, full.names=TRUE)
  
for (fName in q2[8]){
  temp = tryCatch(
    pre1(fName, outDir=outDir),
    warning=function(w) print(paste(w,fName))
  )
  stateTemp = unique(rbind(stateTemp, temp$stateInfo))
  itemTemp = unique(rbind(itemTemp, temp$itemInfo))
}
saveRDS(stateTemp, "./data/state_p13.rds")

2.9 年度資料

上一章節(2.8)的資料是每一年度每一個月的所有測站資料。

此章節主要是要將資料結合成每一年度所有測站的資料。

  • 針對長版資料:

    因為欄位固定,各檔案的欄位均相同。

    將各檔案讀入後,再利用rbind()指令合併即可。

  • 針對寬版資料:

    因各測站的測項並不固定,所以寬版檔案欄位不固定。

    將各檔案讀入後,再利用dplyr::bind_rows()指令合併。

先建立資料夾儲存檔案:

  • S_0165_p_13_year

程式碼如下:

setwd("E:/rumiworkspace")

# 設定輸入資料夾及輸出資料夾
target = "./data/S_0165_P13_G03/"
outDir = "./data/S_0165_p_13_year/"
# 設定長版或寬版
type="W"

fileList = list.files(target, full.names=TRUE)

q1 = unlist(gregexpr(".rds", fileList))

t1 = substr(fileList,q1-8,q1-8)
t2 = substr(fileList,q1-7,q1-4)

yearList = unique(t2[t1==type])

for (year in unique(yearList)){
  fileListTemp = fileList[(t1==type) & (t2==year)]
  
  temp = NULL
  
  for (filename in fileListTemp){
    # 讀入檔案
    # 長版檔案使用
    # 寬版檔案使用
    temp = dplyr::bind_rows(temp, readRDS(filename))
  }
  
  # 寫出
  # 組合檔名
  outFile = paste(outDir, type, year, ".rds", sep="")
  saveRDS(temp, outFile)
}

2.10 補充事項

在進行以上9個章節內容時,可能會發生的錯誤和解決辦法!

(這邊只提較多人出現的情況)

  1. Error in pre1(filename, outDir = outDir) : could not find function “pre1”

Ans:

未創建此函數、未執行此函數、函數名命名與執行寫入函數名不一致。

  1. Error in file(file, “rt”, encoding = fileEncoding) :

    cannot open the connection

Ans:

大多人都是在filename、outDir路徑上出現問題,會去檢查你的資料路徑和名稱是否與你寫入程式一致。

  1. 錯誤: The name of the input file cannot contain the special shell characters: [ <>()|\:&;#?*’] (attempted to copy to a version without those characters ‘Oscar-s-analysis-report.Rmd’ however that file already exists) 停止執行

Ans:

rmarkdown將刪除這些字符以嘗試沒有這些字符的另一個文件。如果我沒記錯的話,我們需要它來正確處理資源。回到儲存該檔案資料夾,將錯誤且多餘的檔案並重新開啟即可。

  1. 以些人在跑大量函數時會遇到,明明其他資料都很順的跑完了,唯獨零星幾個資料跑不出來。

Ans:

我的處理辦法是檢查資料名稱、路徑與寫入程式的是否一致。

若一致,則我認為是資料異常,重新下載資料再執行一次即可。

3.資料庫建置

完成第二章節內容後,我將要把寫好的資料寫入資料庫裡。

3.1 載入RMySQL

library(RMySQL)
## Warning: 套件 'RMySQL' 是用 R 版本 4.2.3 來建造的
## 載入需要的套件:DBI
## Warning: 套件 'DBI' 是用 R 版本 4.2.3 來建造的

3.2 建立資料表

  • 設定連接管道,並將此管道命名為`con`

  • 設定編碼,以避免出現亂碼

  • 設定工作資料庫

    才能將資料放入在我們指定的資料庫。

  • 製作表格名稱

注意:將 - 換成 _底線。

  • 找出適當的SQL語法
可以利用phpmyadmin中的SQL視窗試執行。

CREATE TABLE IF NOT EXISTS 表格名稱(

測站名稱 VARCHAR(3),

測項英文名稱 VARCHAR(15),

測量日期 DATE,

測量小時 TINYINT,

測量值 FLOAT,

測量值old VARCHAR(10),

PRIMARY KEY (測站名稱, 測量日期, 測量小時)

);

再利用套件RMySQL的適當函數執行。

  • 資料插入資料表

  • 關閉連機管道(這個步驟很重要)

否則會創建過多管道,使之後要執行程式會出現下面這錯誤!

如何解決錯誤:

Error in .local(drv, …) : Cannot allocate a new connection: 16 connections already opened

解決辦法:

library(RMySQL)
cons <- dbListConnections(MySQL())
for(con in cons)dbDisconnect(con)

建立資料表的程式碼如下:

userName = "asis"
pswd = "12345678"
hostID = "127.0.0.1"
#source("./myFun/sqlConnectMac.R", local=TRUE)

con = dbConnect(RMySQL::MySQL(), 
                user=userName, 
                password=pswd, 
                host=hostID)

# 設定編碼
strQuery = "SET NAMES UTF8;"
#strQuery = "SET NAMES big5;" # 學校電腦教室需使用 big5
queryResult = dbGetQuery(conn=con, statement=strQuery)

# 設定工作資料庫
# 若無此資料庫,請先建立
strQuery = "USE airp;"
queryResult = dbGetQuery(conn=con, statement=strQuery)

# 讀取資料
setwd("E:/rumiworkspace")
filename = "./data/S_0165_P13_G03/L2014-05.rds"
df = readRDS(filename)

# 製做表格名稱
t1 = unlist(gregexpr(".rds", filename))
tbname = substr(filename, t1-8, t1-1)
# 將 - 換成 _底線
# 資料庫系統中,無法接受『-』為檔案名稱
tbname = gsub("-", "_", tbname)

# 設定查詢語法
strQuery = sprintf("
CREATE TABLE IF NOT EXISTS %s (
  測站名稱 VARCHAR(3),
  測項英文名稱 VARCHAR(15),
  測量日期 DATE,
  測量小時 TINYINT,
  測量值 FLOAT,
  測量值old VARCHAR(10),
  PRIMARY KEY (測站名稱, 測量日期, 測量小時)
);", tbname
)
strQuery
## [1] "\nCREATE TABLE IF NOT EXISTS L2014_05 (\n  測站名稱 VARCHAR(3),\n  測項英文名稱 VARCHAR(15),\n  測量日期 DATE,\n  測量小時 TINYINT,\n  測量值 FLOAT,\n  測量值old VARCHAR(10),\n  PRIMARY KEY (測站名稱, 測量日期, 測量小時)\n);"
# 執行查詢並回收執行結果
queryResult = dbGetQuery(conn=con, statement=strQuery)
# 瀏覽查詢結果
#queryResult

# 關閉連結管道
queryResult = dbDisconnect(conn=con)

3.3 dbWriteTable()

利用RMySQL套件的函數dbWriteTable()插入資料。

函數dbWriteTable()的功用是將R環境中的

  • 資料框(data frame)物件

  • 外部檔案

寫入MySQL資料庫成為資料表,語法如下:

{dbWriteTable(conn, name, value,} field.types=NULL, row.names=TRUE, overwrite=FALSE, append=FALSE, header=TRUE, sep=",", ...)

參數:

  • conn:一個有效的MySQL數據庫連接對象,由dbConnect函數創建。

  • name:要創建的新表的名稱。

  • value:要寫入表中的資料框(data frame)。

  • row.names:一個邏輯值,表示是否將資料框的行名作為MySQL表的一個列。預設值為TRUE,表示將行名作為一個列。

  • field.types:一個可選的參數,用於指定要在MySQL表中創建的各列的數據類型。它是一個命名字符向量,其中鍵是資料框的列名,值是要創建的MySQL列的數據類型。

返回值:

TRUE表示成功將資料框寫入到MySQL表中;

FALSE表示寫入失敗。

3.4 寫成函數

接著我們來接續上面的動作,並寫成函數執行。

tbCreate = function(filename){

  # 載入套件
  library(RMySQL)
  
  # 設定連接管道,並將此管道命名為`con`
  # 注意:帳秘每個人都不一樣
  # host 的設定也有可能不一樣
  userName = "asis"
  pswd = "12345678"
  hostID = "127.0.0.1"
  #source("./myFun/sqlConnectMac.R", local=TRUE)
  
  con = dbConnect(RMySQL::MySQL(), 
                  user=userName, 
                  password=pswd, 
                  host=hostID)
  
  # 設定編碼
  strQuery = "SET NAMES UTF8;"
  #strQuery = "SET NAMES big5;" # 學校電腦教室需使用 big5
  queryResult = dbGetQuery(conn=con, statement=strQuery)
  
  # 設定工作資料庫
  # 若無此資料庫,請先建立
  strQuery = "USE airp;"
  queryResult = dbGetQuery(conn=con, statement=strQuery)
  
  # 讀取資料
  #filename = "./data/aqx_p_13/L2014-01.rds"
  df = readRDS(filename)
  
  # 製做表格名稱
  t1 = unlist(gregexpr(".rds", filename))
  tbname = substr(filename, t1-8, t1-1)
  # 將 - 換成 _底線
  tbname = gsub("-", "_", tbname)

  
  # 資料插入資料表
  bWriteTable(
    conn=con,
    name=tbname,
    value=df,
    row.names=FALSE, 
    overwrite=TRUE,
    append=FALSE
    )
  
  
  # 關閉連結管道
  queryResult = dbDisconnect(conn=con)
}

3.5 大量作業

程式碼與2.8 寫成函數(並大量作業)此一節相似。

註:這邊先以2014年1月作示範,其餘皆有另行下載。

#setwd("E:/rumiworkspace")
#檔案所在目錄(每個人路徑資料夾名稱不同!!!)
#target = "./data/S_0165_P13_G03"
#查出目錄內所有檔案名稱
#q2 = list.files(target, full.names=TRUE)
#找出長版檔案名稱
#type = "L"
#fileList = q2[grepl(type, q2)]
  
# for (fName in fileList[1]){
#    temp = tryCatch(
#    tbCreate(fName),
#     warning=function(w) print(paste(w,fName))
#    )
# }
#}

可能發生的錯誤:

  • Error in .local(conn, statement, …) : could not run statement: Table ‘l2014_06’ already exists

解決辦法:資料來源路徑出現問題,修改路徑

  • cannot open file ‘./myFun/sqlConnectMac.R’: No such file or directoryError in file(filename, “r”, encoding = encoding) : cannot open the connection

解決辦法:#source(“./myFun/sqlConnectMac.R”, local=TRUE)

3.6 補充說明

R中與MySQL數據庫進行連接並操作的優點:

  • 通過與MySQL連接,可以利用R的強大功能對數據進行清洗、轉換、分析和可視化。

  • MySQL是一個廣泛使用的關聯式數據庫管理系統,具有高效的查詢和數據存儲能力。通過將數據輸入到MySQL中,可以利用MySQL的高效能和靈活性進行數據管理、查詢和後續分析。

  • 將數據存儲在MySQL數據庫中,可以實現數據的共享和協作。多個使用者可以通過連接到MySQL數據庫,共同訪問和分析數據,從而促進團隊間的合作。

    例如:陳明輝老師分享給11102-17201統計資料庫班上ASIS資料庫,幫忙學生們都整理好資料,以供我們在[4.研究主題]上好使用。

是否除了RMySQL套件,還有其他一些套件可以實現在R中與MySQL數據庫進行連接並輸入資料的功能或類似之?

  1. RMariaDB:這是另一個在R中與MySQL數據庫交互的套件,提供類似於RMySQL的功能。

    示範:

    library(RMariaDB)

    # 建立與MySQL數據庫的連接

    con <- dbConnect(RMariaDB::MariaDB(), user = “username”, password = “password”,

    dbname = “database_name”, host = “localhost”)

    # 創建一個資料框

    data <- data.frame(ID = c(1, 2, 3), Name = c(“Alice”, “Bob”, “Charlie”))

    # 將資料框寫入到MySQL數據庫中的新表

    dbWriteTable(con, name = “my_table”, value = data)

    # 斷開與MySQL數據庫的連接

    dbDisconnect(con)

  2. RSQLite:這是一個用於在R中與SQLite數據庫進行交互的套件,它提供了類似於RMySQL的函數,如dbWriteTable可以將資料框寫入到SQLite數據庫。

  3. DBI套件:這是一個通用的數據庫接口套件,它提供了一致的API用於在R中與各種數據庫進行交互,包括MySQL、SQLite、PostgreSQL等。通過使用DBI套件,可以使用通用的函數來操作不同的數據庫。

有時在執行 ” library(RMySQL) “ 程式時會出現此提醒

載入需要的套件:DBI

這些套件都提供了類似的函數和功能,可以根據具體的需求和數據庫系統選擇適合的套件。

4.研究主題

4.1 摘要

本研究以環保署空氣品質監測之資料進行試算,討論空氣品質之變化情形。

何謂空氣品質?

空氣品質是指大氣中的空氣成分和污染物的特性和濃度,以及對人類健康和環境的影響程度。

國際上常用的空氣品質指標包括空氣質量指數(Air Quality Index,AQI),就是整合過去空氣污染指標(PSI) 係由每日監測所得懸浮微粒、二氧化硫、一氧化碳、臭氧及二氧化氮等5種主要污染物之濃度值換算而得之空氣污染指標值。細懸浮微粒(PM2.5)指標

AQI值0-50代表空氣品質良好(綠色)51-100普通(黃色)101-150為對敏感族群不良(橘色)151-200為對所有族群不良(紅色)、201-300為非常不良(紫色)301-500有害(褐紅)。從而提供一個簡單明瞭的評估標準,幫助人們了解和應對空氣污染問題。

4.2 研究動機與目的

藉由數據資料,深入了解台灣的空氣品質現狀(所以此主題抓取2022的資料作研究)、污染源以及對人體健康和環境的影響。再以不同方面去探討(所以分成兩大研究主題),通過分析不同地區和城市的空氣質量數據,我們可以評估區域情況並給予建議。

4.3 研究主題一

探討 汙染物的季節變化趨勢 及 地區空氣品質的影響因素分析

4.3.1 子主題一

不同季節下PM2.5和PM10濃度的變化趨勢及分析建議

4.3.2 研究主題說明

根據葉光芃在2019年的10月發表的最新論文指出 ,“粗顆粒帶來的健康效應鮮少受到關注,PM 2.5-10雖無法直達人體內肺泡,卻可滯留在上呼吸道內,會導致發炎和細胞中毒反應,可預見將導致越來越多死亡率,以及呼吸道與心血管疾病的就診率”。

所以此子主題針對pm2.5和pm10這兩測項資料,來看不同季節的變化趨勢,再以不同側站(地區)哪裡是最嚴重(濃度最高)的,深入討論該地區特徵,來得出結論。

4.3.3 方法及結果說明

連接數據庫

userName = "s10170165"
pswd = "s10170165@MCU1234"
hostID = "120.125.72.30"

con = dbConnect(RMySQL::MySQL(), 
                user=userName, 
                password=pswd, 
                host=hostID)
# 設定編碼
strQuery = "SET NAMES UTF8;"
queryResult = dbGetQuery(conn=con, statement=strQuery)
#選擇要用數據庫
dbSendQuery(con, "USE airp15")
## <MySQLResult:878929267,1,1>

獲取日的平均PM2.5和PM10濃度數據(柱狀圖)

SELECT 測量日期, AVG(`PM2.5`), AVG(`PM10`)
FROM W2022_dailyMean
GROUP BY 測量日期
;
Displaying records 1 - 10
測量日期 AVG(PM2.5) AVG(PM10)
2022-01-01 15.46741 27.14263
2022-01-02 16.43269 29.28152
2022-01-03 25.44789 43.46165
2022-01-04 20.65050 34.35375
2022-01-05 24.08404 38.46897
2022-01-06 16.69693 30.79240
2022-01-07 13.95803 26.29919
2022-01-08 15.28793 29.62981
2022-01-09 16.50811 30.11083
2022-01-10 15.38953 27.45594

儲存結果在R中數據框

# 執行SQL查詢
datap1 <- dbGetQuery(con, 
"SELECT 測量日期, AVG(`PM2.5`) 'PM2.5' , AVG(`PM10`) 'PM10'
 FROM W2022_dailyMean
 GROUP BY 測量日期;")

# 打印结果
head(datap1)
##     測量日期    PM2.5     PM10
## 1 2022-01-01 15.46741 27.14263
## 2 2022-01-02 16.43269 29.28152
## 3 2022-01-03 25.44789 43.46165
## 4 2022-01-04 20.65050 34.35375
## 5 2022-01-05 24.08404 38.46897
## 6 2022-01-06 16.69693 30.79240

柱狀圖繪圖

全台日平均PM2.5柱狀圖
library(ggplot2)
## Warning: 套件 'ggplot2' 是用 R 版本 4.2.3 來建造的
# 繪製PM2.5柱狀圖
ggplot(datap1, aes(x = 測量日期, y = PM2.5)) +
  geom_bar(stat = "identity", fill = "blue", width = 0.5) +
  labs(x = "測量日期", y = "平均PM2.5濃度") +
  theme_minimal()

全台日平均PM10柱狀圖
library(ggplot2)

# 繪製PM10柱狀圖
ggplot(datap1, aes(x = 測量日期, y = PM10)) +
  geom_bar(stat = "identity", fill = "green", width = 0.5) +
  labs(x = "測量日期", y = "平均PM10濃度") +
  theme_minimal()

柱狀圖分析

依分析後兩張柱狀圖配合這張全臺細懸浮微粒 (PM2.5) 標準監測數據平均四季變化圖(來源:空氣品質監測網)及PM10歷年月份平均值監測數據變化圖。 可以發現夏秋兩季全臺灣污染物平均濃度較低,而冬季及春季平均濃度較高

獲取季節級別的數據(折線圖)

數據庫中只有日、小時級別的數據而沒有季節級別的數據,我們可以使用SQL查詢和計算來從日級別數據中推導季節級別的數據。以下展示如何計算季節級別的數據:

SELECT
    YEAR(`測量日期`) AS 年份,
    QUARTER(`測量日期`) AS 季度,
    AVG(`PM2.5`) AS 平均PM2_5濃度,
    AVG(`PM10`) AS 平均PM10濃度
FROM
    W2022_dailyMean
GROUP BY
    YEAR(測量日期),
    QUARTER(測量日期)
ORDER BY
    年份,
    季度;
4 records
年份 季度 平均PM2_5濃度 平均PM10濃度
2022 1 18.193016 31.16467
2022 2 11.274614 22.96019
2022 3 9.930115 22.55441
2022 4 13.119943 29.31640

上述代碼使用了YEAR()和QUARTER()函數從日期中提取年份和季度信息,併計算了每個季度的平均PM2.5和PM10濃度。

儲存結果在R中數據框

datap2 <- dbGetQuery(con, 
"SELECT
    YEAR(`測量日期`) AS 年份,
    QUARTER(`測量日期`) AS 季度,
    AVG(`PM2.5`) AS 平均PM2_5濃度,
    AVG(`PM10`) AS 平均PM10濃度
FROM
    W2022_dailyMean
GROUP BY
    YEAR(測量日期),
    QUARTER(測量日期)
ORDER BY
    年份,
    季度;")

折線圖繪圖

全台季節PM2.5變化趨勢圖
library(ggplot2)

# 讀取從SQL查詢得到的季節級別數據
seasonal_data <- datap2 

# 繪製季節級別的PM2.5變化趨勢圖
ggplot(seasonal_data, aes(x = 季度, y = 平均PM2_5濃度)) +
  geom_line(color = "blue") +
  geom_point(color = "blue") +
  labs(x = "2022季度", y = "平均PM2.5濃度") +
  theme_minimal()

全台季節PM10變化趨勢圖
library(ggplot2)

# 讀取從SQL查詢得到的季節級別數據
seasonal_data <- datap2 

# 繪製季節級別的PM10變化趨勢圖
ggplot(seasonal_data, aes(x = 季度, y = 平均PM10濃度)) +
  geom_line(color = "green") +
  geom_point(color = "green") +
  labs(x = "2022季度", y = "平均PM10濃度") +
  theme_minimal()

全台季節PM2.5&PM10變化趨勢圖
ggplot(seasonal_data, aes(x = 季度)) +
  geom_line(aes(y = 平均PM2_5濃度, color = "2.5")) +
  geom_line(aes(y = 平均PM10濃度, color = "PM10")) +
  labs(x = "Date", y = "Value") +
  theme_minimal()

折線圖分析

此折線圖也可以發現夏秋兩季全臺灣污染物平均濃度較低,而冬季及春季平均濃度較高。

而PM10濃度比PM2.5濃度高。

獲取不同站點中PM2.5和PM10濃度最高值(熱力圖)

SELECT
  `測站名稱`,
  MAX(`PM2.5`) AS `最高PM2.5`,
  MAX(`PM10`) AS `最高PM10`
FROM
  `W2022_dailyMean`
GROUP BY
  `測站名稱`
LIMIT 0, 25;
Displaying records 1 - 10
測站名稱 最高PM2.5 最高PM10
三義 33.2917 51.5652
三重 39.2174 57.1739
中壢 43.2609 63.6957
中山 36.1667 52.7917
二林 56.9167 131.4350
仁武 58.9091 80.9167
冬山 32.5000 49.1667
前金 51.9583 81.6667
前鎮 59.1111 95.7222
南投 53.4167 84.2174

儲存結果在R中數據框

datap3 <- dbGetQuery(con, 
"SELECT
  `測站名稱`,
  MAX(`PM2.5`) AS `最高PM2.5`,
  MAX(`PM10`) AS `最高PM10`
FROM
  `W2022_dailyMean`
GROUP BY
  `測站名稱`
LIMIT 0, 25;")

繪圖

PM2.5區域圖

PM2.5排名前五的地區:

嘉義、大寮、前鎮、仁武和善化

# 按照最高PM2.5排序,選擇排名前五的地區
top5 <- head(datap3[order(datap3$最高PM2.5, decreasing = TRUE), ], 5)

# 繪製板塊區域圖
ggplot(top5, aes(x = reorder(測站名稱, -最高PM2.5), y = 1, fill = 最高PM2.5)) +
  geom_tile() +
  labs(x = "測站名稱", y = NULL, fill = "最高PM2.5") +
  scale_fill_gradient(low = "blue", high = "red") +
  theme_minimal() +
  theme(legend.position = "right")

PM10區域圖

PM10排名前五的地區: 二林、大寮、嘉義、前鎮和大里

# 按照最高PM10排序,選擇排名前五的地區
top5 <- head(datap3[order(datap3$最高PM10, decreasing = TRUE), ], 5)

# 繪製板塊區域圖
ggplot(top5, aes(x = reorder(測站名稱, -最高PM10), y = 1, fill = 最高PM10)) +
  geom_tile() +
  labs(x = "測站名稱", y = NULL, fill = "最高PM10") +
  scale_fill_gradient(low = "blue", high = "red") +
  theme_minimal() +
  theme(legend.position = "right")

熱力圖程式碼解釋教學

使用 ggplot() 函數創建繪圖對象,並使用 geom_tile() 函數繪製熱力圖。在 aes() 函數中,將 測站名稱 作為x軸變量,使用1作為y軸變量(此處使用1是為了在圖表中創建一個單獨的行),並使用 最高PM2.5 作為填充顏色。通過 labs() 函數設置坐標軸和圖例的標籤,通過 scale_fill_gradient() 函數設置填充顏色的漸變範圍,最後使用 theme() 函數調整圖表的樣式。

熱力圖分析

根據熱力圖的呈現,可發現大寮、前鎮和嘉義這個測站皆重複於兩張圖上,表示這三個測站的地區在PM2.5和PM10上有些微嚴重,值得去探討其特徵。

大寮和前鎮 是在 高雄的測站

地區解釋(結論)
高雄

大寮和前鎮可能受到工業活動、交通污染和人口密度等因素的影響。

在這些地區,工廠排放、車輛尾氣和大量人口可能導致PM2.5和PM10濃度的升高。

進一步的研究可以關注這些測站所處地區的工業結構、交通流量以及人口分佈情況,以更詳細地了解造成高污染水平的具體原因。

嘉義

嘉義地區顯示出較高的PM2.5和PM10濃度,儘管該地區不屬於高雄地區,但嚴重性也很高。

在乾季(10 月至 4 月)嘉義位於東北季風隻被風處,細懸浮微粒較不容易擴散。

關閉連結通道

queryResult = dbDisconnect(conn=con)

4.3.4 子主題二

各測站空氣汙染物比較

4.3.5 研究主題說明

以上針對完pm2.5和pm10,接著查看各測站汙染物占比,讓數據視覺化,以來看出各地區需特別注意的汙染物預防和趨勢。

這邊的各地區,我們分為北部、中部、雲嘉南、高屏和花東 空品區,再抓取以這五區的地方監測資料的測站綜合值來比較。

4.3.6 方法及結果說明

挑選測站

依據此網站地區中間測站和現有資料測站來挑選

北部:中正、板橋、觀音

中部:大里

雲嘉南:嘉義、崙背、埔里

高屏:鳳山、屏東

花東:花蓮、臺東

連接數據庫

userName = "s10170165"
pswd = "s10170165@MCU1234"
hostID = "120.125.72.30"

con = dbConnect(RMySQL::MySQL(), 
                user=userName, 
                password=pswd, 
                host=hostID)
# 設定編碼
strQuery = "SET NAMES UTF8;"
queryResult = dbGetQuery(conn=con, statement=strQuery)
#選擇要用數據庫
dbSendQuery(con, "USE airp15")
## <MySQLResult:0,2,1>

獲取各地區的地方監測站數據

北部
#北部
SELECT * FROM `W2022_dailyMean` WHERE `測站名稱` IN ('中正','板橋','觀音');
Displaying records 1 - 10
測站名稱 測量日期 AMB_TEMP CH4 CO NMHC NO NO2 NOx O3 PM10 PM2.5 RAINFALL RH SO2 THC WD_HR WIND_DIREC WIND_SPEED WS_HR PH_RAIN RAIN_COND RAIN_INT CO2 Wx.HR Wy.HR Wx.DIRECT Wy.DIRECT
板橋 2022-01-01 17.2375 2.04375 0.398333 0.0875000 1.71250 15.2125 16.9833 29.0333 22.3750 12.00000 0.0000000 79.4583 1.520830 2.13125 120.2500 111.7500 1.04583 0.850000 0.00000 0.00000 0.0000000 NA 0.1142640 0.574366 0.0473125 0.745275
板橋 2022-01-02 17.3833 2.10542 0.430000 0.0979167 2.17917 16.2250 18.4500 24.4667 23.2917 12.45830 0.0000000 76.9167 0.812500 2.20333 96.0417 108.1670 1.27917 1.020830 0.00000 0.00000 0.0000000 NA -0.2328640 0.945135 -0.2918940 1.132690
板橋 2022-01-03 18.8458 2.10000 0.483333 0.0758333 2.52500 16.3583 18.9083 38.6238 34.6250 20.00000 0.0000000 73.2083 0.545833 2.17583 97.8750 102.8750 1.53333 1.141670 0.00000 0.00000 0.0000000 NA -0.0833321 1.071400 -0.1696980 1.349810
板橋 2022-01-04 18.5417 2.01792 0.340417 0.0425000 1.92500 13.5083 15.4792 36.4958 16.1250 9.29167 0.0000000 81.7083 1.837500 2.05792 93.9167 94.0417 2.38750 1.845830 0.00000 0.00000 0.0000000 NA -0.0024818 1.709700 0.0431246 2.178940
板橋 2022-01-05 19.3958 2.11292 0.622083 0.2375000 9.31250 23.8542 33.2083 18.9292 24.4167 14.58330 0.0000000 82.7083 1.470830 2.35042 195.5420 196.1250 1.16667 0.995833 0.00000 0.00000 0.0000000 NA 0.3928380 -0.266902 0.4538090 -0.209908
板橋 2022-01-06 17.1167 2.05316 0.415714 0.1036840 2.12105 17.6105 19.7842 26.9750 19.6522 10.34780 0.3636360 80.7917 0.933333 2.15684 101.0000 103.0870 2.15217 1.608700 1.29167 7.91250 0.3333330 NA -0.3617850 1.410720 -0.4510510 1.846590
板橋 2022-01-07 16.2083 2.01625 0.378696 0.0941667 2.10000 18.2810 20.4143 28.0435 25.3333 11.54170 0.1875000 78.0000 1.065220 2.11042 109.9580 109.0420 1.77500 1.433330 1.09583 2.04167 0.2291670 NA -0.4550400 1.322350 -0.4832120 1.637410
板橋 2022-01-08 16.5042 2.01542 0.316667 0.0420833 1.64167 12.4833 14.1708 32.1333 24.7500 10.25000 0.0000000 74.5833 1.175000 2.05750 103.8330 103.7920 1.79583 1.454170 0.00000 0.00000 0.0000000 NA -0.2655810 1.382720 -0.2864310 1.675800
板橋 2022-01-09 17.8417 2.02417 0.335000 0.0462500 1.25417 12.3750 13.6792 33.2542 21.9167 11.20830 0.0000000 75.6250 1.225000 2.07042 104.3750 112.3750 1.51667 1.233330 0.00000 0.00000 0.0000000 NA -0.1944710 1.147880 -0.3486550 1.360960
板橋 2022-01-10 17.6833 2.01667 0.400000 0.0837500 2.56667 16.9542 19.5708 27.7000 16.9167 9.58333 0.0208333 79.7083 1.000000 2.10042 111.4170 117.5830 1.77083 1.408330 0.23750 2.63750 0.0208333 NA -0.0200464 1.124950 -0.0166594 1.385970
中部
#中部
SELECT * FROM `W2022_dailyMean` WHERE `測站名稱` IN ('大里');
Displaying records 1 - 10
測站名稱 測量日期 AMB_TEMP CH4 CO NMHC NO NO2 NOx O3 PM10 PM2.5 RAINFALL RH SO2 THC WD_HR WIND_DIREC WIND_SPEED WS_HR PH_RAIN RAIN_COND RAIN_INT CO2 Wx.HR Wy.HR Wx.DIRECT Wy.DIRECT
大里 2022-01-01 18.7458 2.23042 0.384583 0.1320830 2.42917 15.5500 18.0458 20.3625 29.0000 18.8333 0.0000000 65.7917 1.62917 2.36250 132.250 178.667 1.77500 1.404170 NA NA NA 444.375 1.114320 0.0006554 1.2305000 -0.0070190
大里 2022-01-02 18.6250 2.28167 0.411250 0.1612500 3.62500 13.7875 17.4625 20.0583 25.8750 17.2083 0.0000000 69.3333 1.61250 2.44292 240.333 253.667 1.60417 1.237500 NA NA NA 450.125 0.658168 -0.1892270 0.8695330 -0.3629620
大里 2022-01-03 19.1167 2.27833 0.559167 0.2208330 8.81250 20.0000 28.8667 23.2000 43.4167 27.2917 0.0000000 70.3333 1.86667 2.49917 265.250 262.000 1.55417 1.225000 NA NA NA 458.417 0.695231 -0.3152430 0.9234620 -0.4058590
大里 2022-01-04 19.4833 2.35083 0.684583 0.4237500 13.53330 26.4750 40.0708 14.8125 51.5000 30.2083 0.0000000 74.9167 2.03750 2.77458 230.458 228.542 1.12917 0.770833 NA NA NA 468.583 0.157299 -0.4042500 0.0710028 -0.4770700
大里 2022-01-05 19.9083 2.34565 0.601739 0.3600000 12.93480 23.6696 36.6435 14.2000 43.0455 26.3913 0.0000000 77.3750 1.76522 2.70565 209.500 206.458 1.33333 0.916667 NA NA NA 464.043 0.511665 -0.2164610 0.6339280 -0.2513420
大里 2022-01-06 18.7208 2.30053 0.416250 0.2357890 4.81250 15.2708 20.1292 19.2583 25.7083 15.7083 0.1666670 76.0833 1.07083 2.53632 223.667 172.750 2.21250 1.700000 NA NA NA 447.917 1.352810 -0.1111270 1.6725000 -0.0492783
大里 2022-01-07 17.5625 2.09429 0.304167 0.1035710 3.37500 12.3208 15.7208 24.1250 14.2083 7.1250 0.1250000 76.7917 1.37500 2.19786 218.125 186.917 2.25833 1.737500 NA NA NA 435.583 1.664270 -0.1720170 1.9106000 -0.2389280
大里 2022-01-08 18.2750 2.13417 0.322083 0.0770833 3.48750 13.0292 16.5667 23.2208 19.6250 10.7500 0.0416667 71.7917 1.65000 2.21125 259.250 291.958 1.92917 1.495830 NA NA NA 436.125 1.389560 -0.3476450 1.7333900 -0.3704840
大里 2022-01-09 18.5583 2.23208 0.358333 0.1120830 5.41250 13.9750 19.4167 24.6542 30.5833 17.6667 0.0000000 74.4583 1.44583 2.34417 230.333 220.958 1.77500 1.441670 NA NA NA 444.792 0.977188 -0.3839520 1.1281800 -0.5592080
大里 2022-01-10 18.5375 2.27542 0.467500 0.1841670 9.29583 17.0875 26.4208 19.6208 35.2083 20.4167 0.0000000 76.3750 1.48750 2.45958 231.792 222.042 1.60000 1.166670 NA NA NA 452.708 0.563592 -0.3770280 0.7574100 -0.5830180
雲嘉南
#雲嘉南
SELECT * FROM `W2022_dailyMean` WHERE `測站名稱` IN ('嘉義','崙背','埔里');
Displaying records 1 - 10
測站名稱 測量日期 AMB_TEMP CH4 CO NMHC NO NO2 NOx O3 PM10 PM2.5 RAINFALL RH SO2 THC WD_HR WIND_DIREC WIND_SPEED WS_HR PH_RAIN RAIN_COND RAIN_INT CO2 Wx.HR Wy.HR Wx.DIRECT Wy.DIRECT
嘉義 2022-01-01 18.7333 2.04500 0.337083 0.121250 1.95417 13.3625 15.3542 25.2167 38.0000 21.8333 0.0000000 68.9583 2.38333 2.16625 45.7083 34.0417 2.63333 2.26250 NA NA NA NA 1.94367 0.9097920 2.14516 1.1544500
嘉義 2022-01-02 18.8125 2.09083 0.361250 0.125000 2.37083 12.2083 14.6208 24.4500 43.3333 25.7083 0.0000000 74.2917 2.58750 2.21583 81.6667 83.0417 1.97917 1.72500 NA NA NA NA 1.57099 0.3733600 1.73073 0.4685170
嘉義 2022-01-03 19.2458 2.13130 0.479565 0.175217 2.83913 18.3391 21.2217 26.7625 57.2273 34.1304 0.0000000 72.5417 3.75217 2.30652 75.2083 74.4167 2.20417 1.92083 NA NA NA NA 1.78946 0.4826990 1.94530 0.5731490
嘉義 2022-01-04 19.7583 2.13833 0.497500 0.227917 3.86250 19.3375 23.2750 25.8095 63.6250 40.6250 0.0000000 76.8750 3.57917 2.36625 143.2920 155.5830 1.60417 1.43333 NA NA NA NA 1.23180 -0.0004508 1.30561 0.0397694
嘉義 2022-01-05 19.8708 2.17042 0.515417 0.252083 4.41250 21.7958 26.2500 15.8417 66.0417 44.2500 0.0000000 83.0833 3.23750 2.42250 93.5417 111.1670 1.81667 1.52917 NA NA NA NA 1.21388 0.4234570 1.47956 0.3906820
嘉義 2022-01-06 19.1625 2.05833 0.378750 0.186667 3.00000 18.1042 21.1583 19.0375 40.7500 23.4167 0.0000000 79.8333 2.78333 2.24500 57.5833 60.0833 2.42917 2.17500 NA NA NA NA 2.07063 0.3788190 2.27950 0.5202170
嘉義 2022-01-07 17.3583 2.01833 0.327083 0.147917 2.77083 16.6333 19.4583 18.8625 26.4167 16.9167 0.0416667 84.4583 2.30833 2.16625 49.5000 78.4167 2.36667 2.07083 NA NA NA NA 1.86465 0.6074310 1.85713 0.5065030
嘉義 2022-01-08 18.3875 2.03917 0.350417 0.131667 2.44167 15.3042 17.7792 23.7875 37.3750 21.2083 0.0000000 75.9583 3.17500 2.17083 91.6667 162.6670 2.04583 1.68750 NA NA NA NA 1.56178 0.3776430 1.87500 0.2891700
嘉義 2022-01-09 18.7458 2.07208 0.341667 0.112500 2.27083 12.7083 15.0042 28.2417 39.4583 24.5417 0.0000000 76.1667 4.16250 2.18458 102.2500 103.7080 2.17917 2.05000 NA NA NA NA 1.93371 0.2691310 1.98095 0.3424650
嘉義 2022-01-10 18.7333 2.07625 0.372083 0.118750 3.10833 13.0000 16.1708 21.9524 41.3913 24.4783 0.0000000 79.7083 2.64167 2.19500 106.0830 127.7500 1.76250 1.57083 NA NA NA NA 1.38046 0.1503190 1.45045 0.2574620
高屏
#高屏
SELECT * FROM `W2022_dailyMean` WHERE `測站名稱` IN ('鳳山','屏東');
Displaying records 1 - 10
測站名稱 測量日期 AMB_TEMP CH4 CO NMHC NO NO2 NOx O3 PM10 PM2.5 RAINFALL RH SO2 THC WD_HR WIND_DIREC WIND_SPEED WS_HR PH_RAIN RAIN_COND RAIN_INT CO2 Wx.HR Wy.HR Wx.DIRECT Wy.DIRECT
屏東 2022-01-01 21.0250 2.15917 0.427917 0.1058330 1.212500 13.6708 14.9292 29.7667 41.6250 25.1250 0.000 59.5417 0.829167 2.26500 254.000 219.542 1.48750 1.22083 NA NA NA NA 0.941790 -0.3061450 1.0974000 -0.2966800
屏東 2022-01-02 21.4958 2.15625 0.372083 0.0820833 0.829167 11.6375 12.5167 33.9708 32.8333 18.4583 0.000 62.5417 1.466670 2.23833 236.875 262.125 1.57083 1.27083 NA NA NA NA 0.553664 -0.3778020 0.5518180 -0.5781180
屏東 2022-01-03 21.8333 2.21125 0.435000 0.1116670 1.704170 14.4875 16.2542 35.5417 45.0000 26.7500 0.000 63.3750 1.433330 2.32292 207.375 176.167 1.47917 1.11667 NA NA NA NA 0.708287 -0.0860183 0.8019740 -0.0812447
屏東 2022-01-04 22.0583 2.21417 0.531667 0.1362500 2.112500 18.1542 20.3167 42.7458 61.2500 37.0833 0.000 66.1250 2.787500 2.35042 167.500 190.917 1.43333 1.16250 NA NA NA NA 0.803308 -0.0967507 0.6945240 -0.1538960
屏東 2022-01-05 23.0042 2.17625 0.424167 0.1108330 0.908333 14.0375 14.9958 33.7750 56.8333 37.2917 0.000 68.5833 1.391670 2.28708 278.375 223.083 1.46250 1.20833 NA NA NA NA 1.054390 -0.3793600 1.1864700 -0.3762350
屏東 2022-01-06 21.5708 2.15875 0.523750 0.1820830 3.645830 19.4292 23.1208 23.2583 64.6250 39.5000 0.000 72.7917 2.816670 2.34083 244.167 275.083 1.45000 1.15417 NA NA NA NA 0.574779 -0.6172840 0.7549260 -0.7416020
屏東 2022-01-07 20.7750 2.09455 0.359545 0.1140910 2.166670 15.8095 18.0143 19.7792 36.0000 22.9048 0.025 74.6667 0.913636 2.20864 288.792 249.208 1.33750 1.11250 NA NA NA NA 0.742434 -0.6652090 0.8584420 -0.6643760
屏東 2022-01-08 21.4667 2.13625 0.438750 0.1425000 2.916670 16.5750 19.5458 35.9714 53.8333 32.6667 0.000 70.8333 2.258330 2.27875 237.125 248.833 1.52083 1.12083 NA NA NA NA 0.156468 -0.6939870 0.0881084 -0.9453870
屏東 2022-01-09 21.3000 2.13375 0.347083 0.0787500 1.404170 10.5000 11.9542 37.1458 47.2500 28.6250 0.000 69.2500 1.675000 2.21250 260.875 282.958 1.66250 1.30000 NA NA NA NA 0.831068 -0.7423670 0.8319360 -1.0849300
屏東 2022-01-10 20.5208 2.12542 0.322917 0.0545833 2.212500 10.2000 12.4583 33.7750 42.5000 25.6250 0.000 72.6250 1.237500 2.18000 290.375 276.750 1.84167 1.43750 NA NA NA NA 1.215780 -0.6249600 1.4704800 -0.7544330
花東
#花東
SELECT * FROM `W2022_dailyMean` WHERE `測站名稱` IN ('花蓮','臺東');
Displaying records 1 - 10
測站名稱 測量日期 AMB_TEMP CH4 CO NMHC NO NO2 NOx O3 PM10 PM2.5 RAINFALL RH SO2 THC WD_HR WIND_DIREC WIND_SPEED WS_HR PH_RAIN RAIN_COND RAIN_INT CO2 Wx.HR Wy.HR Wx.DIRECT Wy.DIRECT
臺東 2022-01-01 19.8333 NA 0.222917 NA 1.437500 4.28750 5.75417 33.6542 11.7083 7.75000 0.0000000 77.3333 1.90833 NA 219.000 232.583 1.99167 1.63333 0.000000 0.000000 0.0000000 NA 1.1244600 -0.1366420 1.352220 -0.1598950
臺東 2022-01-02 20.5167 NA 0.202083 NA 0.658333 3.47500 4.19167 31.2500 13.2917 7.16667 0.0000000 74.0833 1.48750 NA 210.542 238.583 2.05417 1.62083 0.000000 0.000000 0.0000000 NA 1.2885100 -0.2382630 1.605980 -0.3151000
臺東 2022-01-03 20.6417 NA 0.261250 NA 1.352170 3.82174 5.21739 34.9292 17.2609 9.04348 0.0000000 75.2083 1.25833 NA 223.708 238.625 2.20417 1.86250 0.000000 0.000000 0.0000000 NA 1.4070100 -0.3565540 1.677900 -0.3411060
臺東 2022-01-04 20.8125 NA 0.263478 NA 1.117390 4.93478 6.08261 27.7333 12.2174 6.39130 0.0454545 82.0833 1.22174 NA 248.917 237.917 1.40000 1.16250 0.233333 0.358333 0.0208333 NA 0.9126610 -0.3620370 1.070610 -0.3234910
臺東 2022-01-05 22.3417 NA 0.266250 NA 1.141670 5.87083 7.07500 24.1167 13.1250 4.83333 0.0208333 84.8333 1.20000 NA 261.708 236.917 1.37500 0.96250 0.000000 0.000000 0.0000000 NA 0.0660739 -0.3280410 0.169317 -0.3324130
臺東 2022-01-06 20.7292 NA 0.282083 NA 1.537500 5.00417 6.56667 34.8042 12.2500 6.62500 0.0416667 81.0417 1.45417 NA 210.750 238.125 2.12083 1.76250 0.245833 1.070830 0.0208333 NA 1.4885700 -0.1746700 1.775390 -0.2681930
臺東 2022-01-07 19.9792 NA 0.250833 NA 1.216670 4.52917 5.78750 36.2833 10.7917 5.00000 0.0208333 75.6250 1.50417 NA 249.625 248.917 2.15833 1.77083 0.258333 0.587500 0.0208333 NA 1.4142800 -0.3922460 1.706100 -0.4985300
臺東 2022-01-08 19.6833 NA 0.228333 NA 0.912500 4.09167 5.04167 36.7833 10.2917 5.91667 0.0000000 76.6667 1.55833 NA 249.042 236.000 2.19583 1.81250 0.000000 0.000000 0.0000000 NA 1.3981000 -0.4700110 1.650000 -0.4798470
臺東 2022-01-09 20.4500 NA 0.214167 NA 0.850000 3.72083 4.61667 36.1792 14.6667 6.58333 0.0000000 75.9167 1.40000 NA 213.208 226.542 1.88333 1.57500 0.000000 0.000000 0.0000000 NA 1.1995200 -0.2624550 1.460300 -0.3206360
臺東 2022-01-10 21.0458 NA 0.232917 NA 1.695830 3.60833 5.34583 33.9708 15.7500 5.79167 0.0000000 78.7083 1.46250 NA 216.625 230.458 2.12083 1.73333 0.000000 0.000000 0.0000000 NA 1.2490200 0.0381417 1.561820 0.0688152

儲存結果在R中數據框

# 從數據庫中獲取各地區的地方監測站數據
region_data1 <- dbGetQuery(con, "SELECT SUM(`PM2.5`) AS `PM2.5`, SUM(`PM10`) AS `PM10`, SUM(`SO2`) AS `SO2`, SUM(`NO2`) AS `NO2`, SUM(`CO`) AS `CO`, SUM(`O3`) AS `O3` FROM W2022_dailyMean WHERE `測站名稱` IN ('中正','板橋','觀音') GROUP BY `測站名稱`")
region_data1 <- colSums(region_data1)
region_data2 <- dbGetQuery(con, "SELECT  SUM(`PM2.5`) AS `PM2.5`, SUM(`PM10`) AS `PM10`, SUM(`SO2`) AS `SO2`, SUM(`NO2`) AS `NO2`, SUM(`CO`) AS `CO`, SUM(`O3`) AS `O3` FROM W2022_dailyMean WHERE `測站名稱` IN ('大里') GROUP BY `測站名稱`")
region_data2 <- colSums(region_data2)
region_data3 <- dbGetQuery(con, "SELECT  SUM(`PM2.5`) AS `PM2.5`, SUM(`PM10`) AS `PM10`, SUM(`SO2`) AS `SO2`, SUM(`NO2`) AS `NO2`, SUM(`CO`) AS `CO`, SUM(`O3`) AS `O3` FROM W2022_dailyMean WHERE `測站名稱` IN ('嘉義','崙背','埔里') GROUP BY `測站名稱`")
region_data3 <- colSums(region_data3)
region_data4 <- dbGetQuery(con, "SELECT  SUM(`PM2.5`) AS `PM2.5`, SUM(`PM10`) AS `PM10`, SUM(`SO2`) AS `SO2`, SUM(`NO2`) AS `NO2`, SUM(`CO`) AS `CO`, SUM(`O3`) AS `O3` FROM W2022_dailyMean WHERE `測站名稱` IN ('鳳山','屏東') GROUP BY `測站名稱`")
region_data4 <- colSums(region_data4)
region_data5 <- dbGetQuery(con, "SELECT  SUM(`PM2.5`) AS `PM2.5`, SUM(`PM10`) AS `PM10`, SUM(`SO2`) AS `SO2`, SUM(`NO2`) AS `NO2`, SUM(`CO`) AS `CO`, SUM(`O3`) AS `O3` FROM W2022_dailyMean WHERE `測站名稱` IN ('花蓮','臺東') GROUP BY `測站名稱` ;")
region_data5 <- colSums(region_data5)

畫圖

北部
# 可視化各地區的污染物總和
library(ggplot2)
barplot(region_data1,xlab = "北部測站",ylab = "污染物濃度")

中部
# 可視化各地區的污染物總和
library(ggplot2)
barplot(region_data2,xlab = "中部測站",ylab = "污染物濃度")

雲嘉南
# 可視化各地區的污染物總和
library(ggplot2)
barplot(region_data3,xlab = "雲嘉南測站",ylab = "污染物濃度")

高屏
# 可視化各地區的污染物總和
library(ggplot2)
barplot(region_data4,xlab = "高屏測站",ylab = "污染物濃度")

花東
# 可視化各地區的污染物總和
library(ggplot2)
barplot(region_data5,xlab = "花東測站",ylab = "污染物濃度")

結論分析

根據可視化結果,我們可以觀察到所有地區的O3PM10污染物總和相對很高,而PM2.5次高。這可能意味著台灣面臨著較高的臭氧和顆粒物污染問題。

這裡特別提到台北和花東地區測站的PM10也是第二高,但是相對其他地區也明顯的低,那麼我將會在研究主題二有更深入的探討。

高O3濃度:地區可能受到較高的光照和溫暖的氣候影響,這促進了大氣中氮氧化物和揮發性有機化合物的光化反應,產生了更多的臭氧。這種情況通常在夏季和高溫季節更為明顯。

高PM10濃度:地區可能受到交通污染、工業排放和城市化程度高的影響。車輛排放、工業廢氣和揚塵等因素可能導致空氣中顆粒物(PM10)的積累和增加。

這些解釋是基於對污染物總和可視化結果的觀察和常見的空氣污染來源分析。

關閉連結通道

queryResult = dbDisconnect(conn=con)

4.4 研究主題二

都會區(這邊針對台北市)臭氧(O3)每日最高值或惡化情況之分析

4.4.1 子主題一

臭氧濃度與其他空氣汙染物相關性

4.4.2 研究主題說明

首先,我們以環保署所設之空氣品質監測站中的台北市測站為對象,選取2022年1~12月空氣品質測定報告之資料,利用迴歸分析,加以分析及探討。

4.4.3 方法及結果說明

此節為分析架構 與 回歸模型策略

>>取得並整理資料

1. 查詢 W2022_dailyMean 資料表中,有那些是台北市測站(利用排除重複)

程式碼

strQuery = "SELECT DISTINCT 測站名稱 FROM W2022_dailyMean;"

結果
##    測站名稱
## 1      三義
## 2      三重
## 3      中壢
## 4      中山
## 5      二林
## 6      仁武
## 7      冬山
## 8      前金
## 9      前鎮
## 10     南投
## 11     古亭
## 12     善化
## 13     嘉義
## 14     土城
## 15     埔里
## 16     基隆
## 17     士林
## 18     大同
## 19     大園
## 20     大城
## 21     大寮
## 22     大里
## 23     安南
## 24     宜蘭
## 25   富貴角
## 26     小港
## 27     屏東
## 28     崙背
## 29     左營
## 30     平鎮
## 31     彰化
## 32     復興
## 33     忠明
## 34     恆春
## 35     斗六
## 36     新店
## 37     新港
## 38     新營
## 39     新竹
## 40     新莊
## 41     朴子
## 42     松山
## 43     板橋
## 44     林口
## 45     林園
## 46     桃園
## 47     楠梓
## 48     橋頭
## 49     永和
## 50     汐止
## 51     沙鹿
## 52     淡水
## 53     湖口
## 54     潮州
## 55     竹山
## 56     竹東
## 57     線西
## 58     美濃
## 59     臺南
## 60     臺東
## 61     臺西
## 62     花蓮
## 63     苗栗
## 64     菜寮
## 65     萬華
## 66     萬里
## 67     西屯
## 68     觀音
## 69     豐原
## 70     金門
## 71     關山
## 72     陽明
## 73     頭份
## 74     馬公
## 75     馬祖
## 76     鳳山
## 77     麥寮
## 78     龍潭

再根據空氣品質監測網的測站資訊 和 我們排除重複得出的測站資訊,已得到我們之後所針對的測站有:

古亭、士林、大同、陽明、松山、中山、萬華

2. 抓取這些測站的測項資料,並製作成資料框。

程式碼

strQuery = "SELECT `O3`,`CH4`,`CO`,`NMHC`,`SO2`,`PM10`,`PM2.5`,`NO2`,`NO2`,`NO`,`THC`,`CO2`,`CH4` FROM W2022_dailyMean WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')

AND `O3`IS NOT NULL

AND`CH4`IS NOT NULL

AND`CO`IS NOT NULL

AND`NMHC`IS NOT NULL

AND`SO2`IS NOT NULL

AND`PM10`IS NOT NULL

AND`PM2.5`IS NOT NULL

AND`NO2`IS NOT NULL

AND`NO2`IS NOT NULL

AND`NO`IS NOT NULL

AND`THC`IS NOT NULL

AND`CO2`IS NOT NULL

AND`CH4`IS NOT NULL;"

結果
##        O3     CH4       CO     NMHC     SO2    PM10   PM2.5     NO2       NO
## 1 20.3625 2.23042 0.384583 0.132083 1.62917 29.0000 18.8333 15.5500  2.42917
## 2 20.0583 2.28167 0.411250 0.161250 1.61250 25.8750 17.2083 13.7875  3.62500
## 3 23.2000 2.27833 0.559167 0.220833 1.86667 43.4167 27.2917 20.0000  8.81250
## 4 14.8125 2.35083 0.684583 0.423750 2.03750 51.5000 30.2083 26.4750 13.53330
## 5 14.2000 2.34565 0.601739 0.360000 1.76522 43.0455 26.3913 23.6696 12.93480
## 6 19.2583 2.30053 0.416250 0.235789 1.07083 25.7083 15.7083 15.2708  4.81250
##       THC     CO2
## 1 2.36250 444.375
## 2 2.44292 450.125
## 3 2.49917 458.417
## 4 2.77458 468.583
## 5 2.70565 464.043
## 6 2.53632 447.917

3. 檢查是否有 NA 的資料。

程式碼

DATA01[!complete.cases(DATA01),]

結果
##  [1] O3    CH4   CO    NMHC  SO2   PM10  PM2.5 NO2   NO    THC   CO2  
## <0 列> (或零長度的 row.names)

沒有缺失值問題,資料為完整,可以開始進行資料分析。

>>查看資料與繪圖

1. 基本統計量

程式碼

summary(DATA01)

結果
##        O3              CH4              CO              NMHC         
##  Min.   : 3.218   Min.   :1.940   Min.   :0.0987   Min.   :0.007083  
##  1st Qu.:22.021   1st Qu.:2.044   1st Qu.:0.2346   1st Qu.:0.053333  
##  Median :28.654   Median :2.086   Median :0.2946   Median :0.080417  
##  Mean   :28.800   Mean   :2.113   Mean   :0.3236   Mean   :0.106895  
##  3rd Qu.:35.283   3rd Qu.:2.172   3rd Qu.:0.3983   3rd Qu.:0.145217  
##  Max.   :55.383   Max.   :3.374   Max.   :0.7512   Max.   :0.423750  
##       SO2              PM10            PM2.5             NO2        
##  Min.   :0.1958   Min.   : 3.826   Min.   : 0.913   Min.   : 2.608  
##  1st Qu.:0.8917   1st Qu.:13.042   1st Qu.: 6.167   1st Qu.: 5.908  
##  Median :1.2833   Median :20.833   Median :10.250   Median :12.508  
##  Mean   :1.2980   Mean   :23.194   Mean   :12.266   Mean   :12.368  
##  3rd Qu.:1.6042   3rd Qu.:30.625   3rd Qu.:17.000   3rd Qu.:17.279  
##  Max.   :3.6708   Max.   :86.417   Max.   :52.958   Max.   :32.704  
##        NO               THC             CO2       
##  Min.   : 0.1708   Min.   :1.968   Min.   :422.3  
##  1st Qu.: 1.0083   1st Qu.:2.105   1st Qu.:432.0  
##  Median : 1.8583   Median :2.157   Median :437.6  
##  Mean   : 3.0901   Mean   :2.219   Mean   :439.7  
##  3rd Qu.: 3.9542   3rd Qu.:2.314   3rd Qu.:445.9  
##  Max.   :18.7625   Max.   :3.740   Max.   :478.8

2. 變數間相關性

程式碼

cor(DATA01)

結果
##                O3        CH4         CO       NMHC        SO2        PM10
## O3     1.00000000 -0.2653909 -0.3602427 -0.5032774 -0.1836285 -0.03374323
## CH4   -0.26539089  1.0000000  0.6299846  0.7223897  0.4227385  0.52358838
## CO    -0.36024273  0.6299846  1.0000000  0.8306377  0.5256205  0.63252498
## NMHC  -0.50327738  0.7223897  0.8306377  1.0000000  0.5027023  0.59177850
## SO2   -0.18362852  0.4227385  0.5256205  0.5027023  1.0000000  0.53770180
## PM10  -0.03374323  0.5235884  0.6325250  0.5917785  0.5377018  1.00000000
## PM2.5 -0.04750005  0.5756274  0.7108898  0.6252432  0.5474865  0.95764137
## NO2   -0.45871769  0.3582077  0.8172968  0.7171326  0.4916343  0.49951334
## NO    -0.58227987  0.3397323  0.7197607  0.6717651  0.3300462  0.26160190
## THC   -0.38509992  0.9567829  0.7605967  0.8922474  0.4874586  0.59081184
## CO2   -0.47560489  0.5614769  0.7988886  0.7567865  0.4786282  0.59764009
##             PM2.5        NO2         NO        THC        CO2
## O3    -0.04750005 -0.4587177 -0.5822799 -0.3850999 -0.4756049
## CH4    0.57562737  0.3582077  0.3397323  0.9567829  0.5614769
## CO     0.71088980  0.8172968  0.7197607  0.7605967  0.7988886
## NMHC   0.62524317  0.7171326  0.6717651  0.8922474  0.7567865
## SO2    0.54748647  0.4916343  0.3300462  0.4874586  0.4786282
## PM10   0.95764137  0.4995133  0.2616019  0.5908118  0.5976401
## PM2.5  1.00000000  0.5331039  0.3022751  0.6388480  0.6412651
## NO2    0.53310385  1.0000000  0.7640623  0.5352192  0.6976634
## NO     0.30227512  0.7640623  1.0000000  0.5041873  0.6802579
## THC    0.63884797  0.5352192  0.5041873  1.0000000  0.6848578
## CO2    0.64126510  0.6976634  0.6802579  0.6848578  1.0000000

兩兩變數間太相關需要排除!

3. 繪製散佈圖矩陣

程式碼

plot(DATA01)

結果

散佈圖矩陣(Scatterplot Matrix)是一種用於展示多個變數之間關係的圖形展示方式。它以矩陣的形式呈現,每個單元格都是一個散佈圖,顯示了兩個變數之間的關係。在散佈圖矩陣中,對角線上通常是變數的直方圖或密度圖。

散佈圖矩陣的主要目的是讓我們可以快速地視覺化多個變數之間的相互作用和相關性。

4. 變異數膨漲因子

程式碼

下載car套件才可以使用vif函數

install.packages('car') library(car)

DATA01.lm <- lm(Y ~ .,data =DATA01)

vif(DATA01.lm)

結果
## Warning: 套件 'car' 是用 R 版本 4.2.3 來建造的
## 載入需要的套件:carData
## 
## 載入套件:'car'
## 下列物件被遮斷自 'package:dplyr':
## 
##     recode
##          CH4           CO         NMHC          SO2         PM10        PM2.5 
## 5.667093e+04 8.000079e+00 2.358244e+04 1.602005e+00 1.390121e+01 1.745729e+01 
##          NO2           NO          THC          CO2 
## 4.690117e+00 3.681664e+00 1.330272e+05 3.440683e+00

*判斷多元線性迴歸模型的自變數之間是否獨立,VIF值越小越好,若VIF值>10,表示自變數存在共線性,則應刪除該自變數

所以我們應該刪除 CH4、NMHC、PM10、PM2.5、THC 這些變數(測項)

DATA01n <- DATA01[,-c(2,4,6,7,10)]
head(DATA01n)
##        O3       CO     SO2     NO2       NO     CO2
## 1 20.3625 0.384583 1.62917 15.5500  2.42917 444.375
## 2 20.0583 0.411250 1.61250 13.7875  3.62500 450.125
## 3 23.2000 0.559167 1.86667 20.0000  8.81250 458.417
## 4 14.8125 0.684583 2.03750 26.4750 13.53330 468.583
## 5 14.2000 0.601739 1.76522 23.6696 12.93480 464.043
## 6 19.2583 0.416250 1.07083 15.2708  4.81250 447.917

>>建立迴歸模型

1. 建立迴歸模型

DATA01n.lm <- lm(O3 ~ CO+SO2+NO2+NO+CO2 ,data = DATA01n)

2. 查看模型配適結果(含參數檢定)

程式碼

summary(DATA01n.lm)

結果
## 
## Call:
## lm(formula = O3 ~ CO + SO2 + NO2 + NO + CO2, data = DATA01n)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -21.4604  -4.8924  -0.0238   4.0438  22.4849 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 164.77882   25.42172   6.482 2.58e-10 ***
## CO           33.98468    6.43558   5.281 2.08e-07 ***
## SO2           0.26734    0.81873   0.327   0.7442    
## NO2          -0.25385    0.11223  -2.262   0.0242 *  
## NO           -1.60102    0.19246  -8.319 1.29e-15 ***
## CO2          -0.31667    0.06066  -5.220 2.83e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.54 on 415 degrees of freedom
## Multiple R-squared:  0.3943, Adjusted R-squared:  0.387 
## F-statistic: 54.03 on 5 and 415 DF,  p-value: < 2.2e-16

解釋: Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ’ ’ 1

是標註預測變量的顯著性”***“最顯著,然後依次減小

(由於變數顯著量過少,所以判斷為有些變數影響了其他變數的顯著性)

3. 減少變數

因此使用Forward 和 Backward(AIC準則)

Forward(前向選擇)

是一種逐步添加自變量的方法。它從一個空模型開始,然後在每一步選擇與目標變量之間具有最高相關性的一個自變量,直到滿足某個停止準則(如預先設定的顯著性水平)或沒有其他可添加的自變量為止。這種方法逐步增加自變量,每一步都根據模型的性能指標(例如,R方或AIC)來評估加入自變量後模型的改進程度。

程式碼

null = lm(Y ~ 1, data = DATA01n)

full = lm(Y ~ ., data = DATA01n)

forward.lm = step(null,scope = list(lower=null, upper=full),direction = "forward")

summary(forward.lm)

結果
## Start:  AIC=1908.03
## O3 ~ 1
## 
##        Df Sum of Sq   RSS    AIC
## + NO    1   13205.8 25744 1735.7
## + CO2   1    8810.4 30139 1802.1
## + NO2   1    8195.8 30754 1810.6
## + CO    1    5054.7 33895 1851.5
## + SO2   1    1313.4 37636 1895.6
## <none>              38950 1908.0
## 
## Step:  AIC=1735.7
## O3 ~ NO
## 
##        Df Sum of Sq   RSS    AIC
## + CO2   1    458.26 25286 1730.1
## + CO    1    279.99 25464 1733.1
## <none>              25744 1735.7
## + NO2   1     17.87 25726 1737.4
## + SO2   1      3.20 25741 1737.7
## 
## Step:  AIC=1730.14
## O3 ~ NO + CO2
## 
##        Df Sum of Sq   RSS    AIC
## + CO    1   1402.90 23883 1708.1
## <none>              25286 1730.1
## + SO2   1    107.64 25178 1730.3
## + NO2   1     17.05 25268 1731.9
## 
## Step:  AIC=1708.11
## O3 ~ NO + CO2 + CO
## 
##        Df Sum of Sq   RSS    AIC
## + NO2   1    285.15 23597 1705.1
## <none>              23883 1708.1
## + SO2   1      0.39 23882 1710.1
## 
## Step:  AIC=1705.06
## O3 ~ NO + CO2 + CO + NO2
## 
##        Df Sum of Sq   RSS    AIC
## <none>              23597 1705.1
## + SO2   1     6.061 23591 1707.0
## 
## Call:
## lm(formula = O3 ~ NO + CO2 + CO + NO2, data = DATA01n)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -21.5643  -4.8498  -0.0709   4.0718  22.5923 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 163.75095   25.19897   6.498 2.33e-10 ***
## NO           -1.61242    0.18906  -8.528 2.78e-16 ***
## CO2          -0.31392    0.06001  -5.231 2.67e-07 ***
## CO           34.35265    6.32934   5.428 9.73e-08 ***
## NO2          -0.24724    0.11027  -2.242   0.0255 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.532 on 416 degrees of freedom
## Multiple R-squared:  0.3942, Adjusted R-squared:  0.3883 
## F-statistic: 67.66 on 4 and 416 DF,  p-value: < 2.2e-16

Backward(後向選擇)

是一種逐步剔除自變量的方法。它從包含所有自變量的完整模型開始,然後在每一步中,剔除對模型性能貢獻較小的一個自變量,直到滿足某個停止準則(如預先設定的顯著性水平)或沒有其他可剔除的自變量為止。這種方法逐步減少自變量,每一步都根據模型的性能指標來評估去除自變量後模型的改進程度。

程式碼

full = lm(O3 ~ ., data = DATA01n)

backward.lm = step(full, scope = list(upper=full),direction = "backward")

summary(backward.lm)

結果
## Start:  AIC=1706.95
## O3 ~ CO + SO2 + NO2 + NO + CO2
## 
##        Df Sum of Sq   RSS    AIC
## - SO2   1       6.1 23597 1705.1
## <none>              23591 1707.0
## - NO2   1     290.8 23882 1710.1
## - CO2   1    1549.2 25141 1731.7
## - CO    1    1585.2 25177 1732.3
## - NO    1    3933.9 27525 1769.9
## 
## Step:  AIC=1705.06
## O3 ~ CO + NO2 + NO + CO2
## 
##        Df Sum of Sq   RSS    AIC
## <none>              23597 1705.1
## - NO2   1     285.1 23883 1708.1
## - CO2   1    1552.3 25150 1729.9
## - CO    1    1671.0 25268 1731.9
## - NO    1    4125.9 27723 1770.9
## 
## Call:
## lm(formula = O3 ~ CO + NO2 + NO + CO2, data = DATA01n)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -21.5643  -4.8498  -0.0709   4.0718  22.5923 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 163.75095   25.19897   6.498 2.33e-10 ***
## CO           34.35265    6.32934   5.428 9.73e-08 ***
## NO2          -0.24724    0.11027  -2.242   0.0255 *  
## NO           -1.61242    0.18906  -8.528 2.78e-16 ***
## CO2          -0.31392    0.06001  -5.231 2.67e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.532 on 416 degrees of freedom
## Multiple R-squared:  0.3942, Adjusted R-squared:  0.3883 
## F-statistic: 67.66 on 4 and 416 DF,  p-value: < 2.2e-16

4.精煉模型

根據原始模型、Forward 和 Backward(AIC準則的結果,可以一一推出結論。

  • 先刪除不必要的

如:NO2(顯著不高移除)、NO(下面解釋)

  • 新增已刪變數(是否可提高Multiple R-squared and Adjusted R-squared)

如:PM10(Multiple R-squared and Adjusted R-squared 皆有提升)

小問題1 :這樣會顯得NO顯著不高,但刪除NO變向會使Multiple R-squared and Adjusted R-squared 降低,這樣解釋模型的能力會太低,所以還是保留NO!

小問題2 :當NO留下則會使CO顯著性降低(>0.05),除我們刪除CO

最後可得出

DATA01x.lm <- lm(O3 ~ CO2+NO+PM10 ,data = DATA01)
summary(DATA01x.lm)
## 
## Call:
## lm(formula = O3 ~ CO2 + NO + PM10, data = DATA01)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -24.4916  -4.5140   0.1475   4.0131  20.8833 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 181.01648   24.94121   7.258 1.94e-12 ***
## CO2          -0.34935    0.05858  -5.963 5.28e-09 ***
## NO           -1.20827    0.16030  -7.538 3.01e-13 ***
## PM10          0.22084    0.03588   6.156 1.76e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.456 on 417 degrees of freedom
## Multiple R-squared:  0.4049, Adjusted R-squared:  0.4006 
## F-statistic: 94.57 on 3 and 417 DF,  p-value: < 2.2e-16

5. 模型定型

結論

得出:O3 ~ CO2+NO+PM10

DATA01x.lm <- lm(O3 ~ CO2+NO+PM10 ,data = DATA01)

由此得知O3與 CO2、NO、PM10 這些汙染物有相關性。

圖形
# 加載所需的套件
library(ggplot2)

# 擬合回歸模型
model <- lm(O3 ~ CO2+NO+PM10 ,data = DATA01)
DATA01$predicted <- predict(model)

# 創建回歸模型圖
ggplot(DATA01, aes(x = predicted, y = O3)) +
  geom_point() +                        # 繪制散點圖
  geom_smooth(method = "lm", se = FALSE) # 添加回歸線,se = FALSE禁用置信區間
## `geom_smooth()` using formula = 'y ~ x'

上面的模型中,我們使用數據集DATA01,其中CO2、NO和PM10是三個自變量,O3是因變量。然後,我們使用lm()函數擬合了一個多元線性回歸模型model,其中O3 ~ CO2+NO+PM10表示使用CO2、NO和PM10預測O3。

接下來,我們使用predict()函數計算預測值,並將預測值添加到數據集中的新列predicted。最後,使用ggplot2創建散點圖,並使用geom_smooth()函數添加擬合線。

6. 模型診斷與檢定

畫出模型診斷用的圖

程式碼

需先載入ggplot2、ggfortify才可以使用autoplot函數

library(ggplot2)

library(ggfortify)

autoplot(DATA01x.lm)

結果
## Warning: 套件 'ggfortify' 是用 R 版本 4.2.3 來建造的

Residuals vs. Fitted

表現還行的殘差圖

Normal QQ

符合常態分配

Scale-Location

這是比例位置圖,它沿 x 軸顯示回歸模型的擬合值,沿 y 軸顯示標準化殘差的平方根

然而此圖殘差的分佈在所有擬合值處大致相等

Residuals vs Leverage

這是殘差與槓桿圖,是一種診斷圖,它能夠識別回歸模型中有影響的觀察結果。

我們的Residuals vs Leverage圖中,表示我們的回歸模型中沒有任何影響點

>>模型檢定能力

1. 殘差常態性檢定

程式碼

shapiro.test(DATA01x.lm$residual)

結果
## 
##  Shapiro-Wilk normality test
## 
## data:  DATA01x.lm$residual
## W = 0.99349, p-value = 0.06664

常態性檢定的 p-value 以一般 95% 的信賴水準來說,是不拒絕虛無假設的,也就是說殘差值的常態分配假設是合理的。

2. 殘差獨立性檢定

程式碼

durbinWatsonTest(DATA01x.lm)

結果
##  lag Autocorrelation D-W Statistic p-value
##    1       0.6102951     0.7757695       0
##  Alternative hypothesis: rho != 0

從輸出中我們可以看到檢驗統計量為 0.7757695,對應的 p 值為 0。由於此 p 值為 0,我們可以拒絕零假設並得出結論,此回歸模型中的殘差是完全正自相關的。

3. 殘差變異數同質性檢定

程式碼

ncvTest(DATA01x.lm)

結果
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 2.138784, Df = 1, p = 0.14362

這裡發現殘差變異數同質性檢定,以一般 95% (0.05)的信賴水準來說,p = 0.14362 是不拒絕虛無假設的,也就是說殘差的變異數有符合同質性的假設。

>>結論

最後我們的迴歸模型為O3 ~ -0.34935(CO2)-1.20827(NO)+0.22084(PM10)+181.01648

由上述模型可得知當PM10下降時 O3也會隨之下降CO2、NO上升時 O3下降,所以我們可以根據此模型來算出地區環境汙染物的平衡點,是對於該區域是最佳的!

4.4.4 子主題二

臭氧和高影響污染物的資料視覺化

4.4.5 研究主題說明

可以發現PM10的測值高時也相對造成臭氧測值的偏高,所以此主題以圖表展現,讓數據呈現減少那些汙染物可以抑制臭氧。

4.4.6 方法及結果說明

我們使用R和RMySQL套件來製作圖表,並從MySQL資料庫中擷取相關數據。以下是程式碼和結果,展示如何使用RMySQL套件連接到MySQL資料庫並製作圖表。

>>抓取資料

先從資料庫抓取所需資訊(O3&PM10),並命名此資料表為DATA02

library(RMySQL)

# 設定資料庫連線資訊
con <- dbConnect(RMySQL::MySQL(), 
                 user = "s10170165", 
                 password = "s10170165@MCU1234", 
                 host = "120.125.72.30")

# 設定編碼
strQuery = "SET NAMES UTF8;"
queryResult = dbGetQuery(conn=con, statement=strQuery)

# 設定工作資料庫
# 若無此資料庫,請先建立
strQuery = "USE airp15;"
queryResult = dbGetQuery(conn=con, statement=strQuery)

# 執行SQL查詢並擷取數據
strQuery <- "SELECT `PM10`, `O3` 
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL;"
DATA02 <- dbGetQuery(conn=con, statement=strQuery)

# 關閉資料庫連線
queryResult = dbDisconnect(conn=con)

#檢查資料
head(DATA02)
##      PM10      O3
## 1 16.8750 29.6208
## 2 16.7083 33.2708
## 3 35.2083 39.8417
## 4 20.3913 28.2750
## 5 22.6087 20.9958
## 6 15.0417 31.7667

>>資料整理

這邊要以季、月做一次整理,否則在做一些圖時,以日的資料太龐大,畫出來的圖不好分析或是圖表呈現毫無意義!

抓取四季資料並建置好表。

library(RMySQL)

# 設定資料庫連線資訊
con <- dbConnect(RMySQL::MySQL(), 
                 user = "s10170165", 
                 password = "s10170165@MCU1234", 
                 host = "120.125.72.30")

# 設定編碼
strQuery = "SET NAMES UTF8;"
queryResult = dbGetQuery(conn=con, statement=strQuery)

# 設定工作資料庫
# 若無此資料庫,請先建立
strQuery = "USE airp15;"
queryResult = dbGetQuery(conn=con, statement=strQuery)

# 執行SQL查詢並擷取數據
strQuery <- "SELECT `PM10`, `O3` 
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)<=3
;"
DATA_1 <- dbGetQuery(conn=con, statement=strQuery)
strQuery <- "SELECT `PM10`, `O3` 
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)<=6
AND MONTH(測量日期)>=4
;"
DATA_2 <- dbGetQuery(conn=con, statement=strQuery)
strQuery <- "SELECT `PM10`, `O3` 
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)<=9
AND MONTH(測量日期)>=7
;"
DATA_3 <- dbGetQuery(conn=con, statement=strQuery)
strQuery <- "SELECT `PM10`, `O3` 
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)<=12
AND MONTH(測量日期)>=10
;"
DATA_4 <- dbGetQuery(conn=con, statement=strQuery)


# 關閉資料庫連線
queryResult = dbDisconnect(conn=con)

接著整理季份資料,變成占比(季/年)

b <- colSums(DATA02)

a1 <- colSums(DATA_1)
Quarter1 <- a1 / b

a2 <- colSums(DATA_2)
Quarter2 <- a2 / b

a3 <- colSums(DATA_3)
Quarter3 <- a3 / b

a4 <- colSums(DATA_4)
Quarter4 <- a4 / b

d <- data.frame(Quarter1, Quarter2, Quarter3,  Quarter4)
d
##       Quarter1  Quarter2  Quarter3  Quarter4
## PM10 0.2819162 0.2187419 0.2167358 0.2826061
## O3   0.2593148 0.2324895 0.2398352 0.2683605
# 转置数据
e <- data.frame(t(d))
e
##               PM10        O3
## Quarter1 0.2819162 0.2593148
## Quarter2 0.2187419 0.2324895
## Quarter3 0.2167358 0.2398352
## Quarter4 0.2826061 0.2683605
# 增加季度
Quarter <- 1:4
f <- cbind(e, Quarter)
f
##               PM10        O3 Quarter
## Quarter1 0.2819162 0.2593148       1
## Quarter2 0.2187419 0.2324895       2
## Quarter3 0.2167358 0.2398352       3
## Quarter4 0.2826061 0.2683605       4

抓取單月份資料

library(RMySQL)

# 設定資料庫連線資訊
con <- dbConnect(RMySQL::MySQL(), 
                 user = "s10170165", 
                 password = "s10170165@MCU1234", 
                 host = "120.125.72.30")

# 設定編碼
strQuery = "SET NAMES UTF8;"
queryResult = dbGetQuery(conn=con, statement=strQuery)

# 設定工作資料庫
# 若無此資料庫,請先建立
strQuery = "USE airp15;"
queryResult = dbGetQuery(conn=con, statement=strQuery)

# 執行SQL查詢並擷取數據
strQuery <- "SELECT `PM10`, `O3` 
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=1
;"
DATA_101 <- dbGetQuery(conn=con, statement=strQuery)

strQuery <- "SELECT `PM10`, `O3` 
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=2
;"
DATA_102 <- dbGetQuery(conn=con, statement=strQuery)

strQuery <- "SELECT `PM10`, `O3` 
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=3
;"
DATA_103 <- dbGetQuery(conn=con, statement=strQuery)

strQuery <- "SELECT `PM10`, `O3` 
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=4
;"
DATA_104 <- dbGetQuery(conn=con, statement=strQuery)

strQuery <- "SELECT `PM10`, `O3` 
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=5
;"
DATA_105 <- dbGetQuery(conn=con, statement=strQuery)

strQuery <- "SELECT `PM10`, `O3` 
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=6
;"
DATA_106 <- dbGetQuery(conn=con, statement=strQuery)

strQuery <- "SELECT `PM10`, `O3` 
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=7
;"
DATA_107 <- dbGetQuery(conn=con, statement=strQuery)

strQuery <- "SELECT `PM10`, `O3` 
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=8
;"
DATA_108 <- dbGetQuery(conn=con, statement=strQuery)

strQuery <- "SELECT `PM10`, `O3` 
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=9
;"
DATA_109 <- dbGetQuery(conn=con, statement=strQuery)

strQuery <- "SELECT `PM10`, `O3` 
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=10
;"
DATA_110 <- dbGetQuery(conn=con, statement=strQuery)

strQuery <- "SELECT `PM10`, `O3` 
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=11
;"
DATA_111 <- dbGetQuery(conn=con, statement=strQuery)

strQuery <- "SELECT `PM10`, `O3` 
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=12
;"
DATA_112 <- dbGetQuery(conn=con, statement=strQuery)


# 關閉資料庫連線
queryResult = dbDisconnect(conn=con)

接著整理月份資料

一月<-colSums(DATA_101)
二月<-colSums(DATA_102)
三月<-colSums(DATA_103)
四月<-colSums(DATA_104)
五月<-colSums(DATA_105)
六月<-colSums(DATA_106)
七月<-colSums(DATA_107)
八月<-colSums(DATA_108)
九月<-colSums(DATA_109)
十月<-colSums(DATA_110)
十一月<-colSums(DATA_111)
十二月<-colSums(DATA_112)
m = data.frame(一月,二月,三月,四月,五月,六月,七月,八月,九月,十月,十一月,十二月)
m
##          一月     二月     三月     四月     五月     六月     七月     八月
## PM10 75556.85 53421.99 80852.89 74094.79 48673.03 40042.91 44020.20 44567.90
## O3   69644.20 66384.30 76538.31 78732.08 70773.79 41071.56 50187.53 56831.98
##          九月     十月   十一月   十二月
## PM10 72729.45 71454.73 72876.52 66013.94
## O3   89579.37 81402.93 68808.98 69769.89
#轉置資料
m<-data.frame(t(m))
m
##            PM10       O3
## 一月   75556.85 69644.20
## 二月   53421.99 66384.30
## 三月   80852.89 76538.31
## 四月   74094.79 78732.08
## 五月   48673.03 70773.79
## 六月   40042.91 41071.56
## 七月   44020.20 50187.53
## 八月   44567.90 56831.98
## 九月   72729.45 89579.37
## 十月   71454.73 81402.93
## 十一月 72876.52 68808.98
## 十二月 66013.94 69769.89
#增加季度
<-1:12
m <- cbind(m,月)
m
##            PM10       O3 月
## 一月   75556.85 69644.20  1
## 二月   53421.99 66384.30  2
## 三月   80852.89 76538.31  3
## 四月   74094.79 78732.08  4
## 五月   48673.03 70773.79  5
## 六月   40042.91 41071.56  6
## 七月   44020.20 50187.53  7
## 八月   44567.90 56831.98  8
## 九月   72729.45 89579.37  9
## 十月   71454.73 81402.93 10
## 十一月 72876.52 68808.98 11
## 十二月 66013.94 69769.89 12

>> 製圖

散布圖

程式碼
# 載入ggplot2套件
library(ggplot2)

# 製作散點圖
ggplot(DATA02, aes(x = PM10, y = O3)) +
  geom_point() +
  xlab("PM10") +
  ylab("O3")
結果

散點圖展示了PM10值和臭氧值之間的關係。在這個圖表中,每個散點代表一個數據點,其中x軸表示PM10值,y軸表示臭氧值。

折線圖

程式碼
# 載入ggplot2套件
library(ggplot2)

# 繪製折線圖
ggplot() +
  geom_line(data=f, aes(x = Quarter, y = PM10,color = "PM10"),size=1) +
  geom_point(data=f, aes(x = Quarter, y = PM10,color = "PM10"),size=3)+
  
  geom_line(data=f, aes(x = Quarter, y = O3,color = "O3"),size=1) +
  geom_point(data=f, aes(x = Quarter, y = O3,color = "O3"),size=3)+
  xlab("季度")+ylab("PM10 & O3")
結果
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

上述折線圖展示了PM10和臭氧值隨季節變化的趨勢。這種圖表通常用於觀察兩個變數之間的關係,並了解它們隨著時間的變化而產生的趨勢。

在折線圖中,x軸表示季度,而y軸表示數值變量,例如PM10值和臭氧值。每條折線代表一個變量的變化。如果有多個變量,它們可以使用不同的顏色或線型來區分,像這邊我們O3為紅、PM10為藍。

通過觀察折線的趨勢,我們可以得出一些初步的解釋:

  1. 趨勢方向:觀察折線的整體趨勢,我們可以發現這兩汙染物在第二、三季是處於低處,而在第一、四季處於高點

  2. 趨勢變化率:觀察折線的斜率變化,可以了解PM10和臭氧值變化的速率。陡峭的斜率表示變化快速,而平緩的斜率表示變化緩慢。

  3. 相對位置:觀察兩條折線之間的相對位置,可以了解PM10和臭氧值之間的時間延遲關係。如果PM10的變化在臭氧值之前發生,則可以推斷PM10對臭氧的影響有一定的時間滯後。

需要注意的是,台灣PM10、O3濃度之變化除了兩者互相影響外,之變化也易受於季節影響,主要約集中於10月至2月

柱狀圖

各別汙染物
# 載入ggplot2套件
library(ggplot2)

# 製作各別柱狀圖
ggplot(m, aes(x = 月, y = PM10)) +
  geom_bar(stat = "identity", position = "dodge") +
  xlab("月") +
  ylab("PM10")

ggplot(m, aes(x = 月, y = O3)) +
  geom_bar(stat = "identity", position = "dodge") +
  xlab("月") +
  ylab("O3")

在柱狀圖中,x軸表示月,而y軸個各別表示臭氧值、PM10。

可以觀察不同月份的臭氧值、PM10柱子的變化趨勢。

O3 vs PM10
# 載入ggplot2套件
library(ggplot2)

# 繪製柱狀圖
ggplot(DATA02, aes(x = PM10, y = O3)) +
  geom_bar(stat = "identity", position = "dodge") +
  xlab("PM10") +
  ylab("O3")

可以推斷PM10對臭氧的影響是正向還是負向的。如果隨著PM10水平的增加,臭氧值的柱子也增加,則表示PM10對臭氧具有促進作用;反之,如果臭氧值的柱子隨著PM10水平的增加而減少,則表示PM10對臭氧具有抑制作用。

結論

根據研究結果,減少PM10的排放可以有效抑制臭氧(O3)的生成,這在都會區(北部)的環境管理中扮演著重要的角色。

研究選取了台北市環保署設立的七個空氣品質監測站作為研究對象,分析了2022年1至12月的空氣品質測定報告資料,並使用迴歸分析方法進行了相關性分析和探討。

研究結果顯示,在影響臭氧生成的因子中,CO2、NO和PM10對臭氧濃度的影響最大。透過迴歸分析,我們得知O3濃度的變化可以使用CO2、NO和PM10來預測,這三個因子解釋了48%的變異量(R^2 = 0.4899)。

從這些結果中,我們可以觀察到當PM10的濃度較高時,臭氧的濃度也相對較高。因此,如何有效控制這些污染物的生成量,特別是針對對臭氧生成影響最大的污染物的減排,被列為抑制臭氧生成的首要目標

總結而言,這項研究的結果強調了北部減少PM10排放對於抑制臭氧生成的重要性。進一步的研究和環境管理措施可以針對高影響臭氧生成的污染物,尤其是PM10,以減少其排放量,從而有效改善都會區的空氣品質。